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The  dynamics  of  ocean  fronts  is  an  important  aspect  of  ocean 
theory  which  has  been  overlooked  until  recently.  The  excellent 
observational  work  of  Gunner  Roden  has  provided  us  with  very  useful 
data  on  which  to  build  models.  The  present  work  is  a  rather  simple 
model  of  the  Pacific  Ocean  Sub-Arctic  Front  forced  by  atmospheric 
wind  events. 

Most  of  the  results  of  this  work  advance  computational  oceanography 
rather  than  shed  new  insight  on  ocean  front  dynamics.  We  have  learned 
how  to  formulate  an  isopycnal  coordinate  system  for  this  problem; 
how  to  design  the  open  boundary  conditions  and  how  to  specify  the 
initial  conditions.  An  old  problem  of  instability  in  layered  systems 
is  resolved  by  recognizing  that  the  eigenvectors  of  the  layered 
system  may  introduce  an  anomalous  fast  mode  if  the  parameters  are 
not  properly  chosen. 

It  is  shown  analytically  that  instability  due  to  horizontal  shear 
is  only  possible  on  very  small  scales  and  should  not  be  active  in  the 
sub-arctic  front.  In  addition  we  see  that  the  wind  systems  only 
influence  the  upper  ocean  and  do  not  penetrate  to  very  deep  layers 
except  with  extreme  forcing.  This  is  consistent  with  the  observations 
of  Roden. 

This  research  is  only  a  beginning  in  our  effort  to  model  upper 
ocean  fronts.  Considerable  more  work  is  required. 

James  J.  O'Brien 

Director 

Mesoscale  Air-Sea  Interaction  Group 
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Theoretical  studies  and  numerical  experiments  of  the  subarctic  front 
are  conducted.  The  theoretical  side  of  this  research  encompasses  the 
determination  of  a  stability  criterion  and  the  role  friction  plays  as  a 
dissipative  mechanism.  In  the  first  case  study,  it  is  established  that 
if  the  difference  between  the  phase  speed  of  the  perturbation  and  the 
mean  long- front  geostrophic  flow  is  greater  than  a  certain  value,  the 
subarctic  front  is  stable  everywhere.  In  the  other  case  study,  it  is 
proved  that  the  time  scale  of  dissipation  is  sufficiently  large  to 
neglect  horizontal  friction.  While  the  effect  of  the  westerly  winds 
causes  a  southward  shift  of  the  front,  the  effect  of  a  negative  curl 
shows  convergence  at  the  front.  This  convergence  strongly  favors 
froncogenesis. 

It  is  established  that  the  phase  speed  of  the  first  baroclinic  mode 
is  larger  in  a  numerical  model  where  a  constant  density  jump  between 
layers  is.  employed.  Thus,  a  smaller  time  scale  is  required  to  avoid  the 
numerical  instability  caused  by  the  violation  of  the  CFL  criterion.  A 
numerical  four-layer  hydrodynamics!  model  is  implemented  to  investigate 
the  role  the  atmosphere  exerts  on  the  dynamics  of  the  subarctic  front. 

The  salient  features  at  the  wake  of  a  strong  cyclone  can  be 
summarized  as  follows:  (i)  upwelling  of  20  meters  of  the  four  interfaces 
occurs;  (ii)  the  water  remains  upwellad  for  several  weeks;  (iii)  the 
e-folding  time  scale  of  the  vertical  oscillations  induced  by  the  storm's 
passage  is  of  the  order  of  10  aays;  (iv)  a  region  of  relative  maximum 
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convergence  is  observed  in  the  lower  layers  below  a  region  of  relative 
maximum  upwelling  in  the  upper  layer;  and  (v)  the  long  and  cross-front 
velocity  fields  are  in  quadrature  In  time,  with  the  cross  front  velocity 
lagging  the  long  front  velocity. 

With  the  exception  of  the  cyclone's  case,  it  is  established  that  the 
effects  of  the  atmospheric  wind  forcing  are  confined  to  the  upper  layer. 
These  results  are  supported  by  oceanic  observations.  Finally,  it  is  also 
proved  that  the  response  of  the  front  to  the  atmospheric  momentum 
transfer  is  independent  of  the  initial  form  of  the  subarctic  front. 
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1.  INTRODUCTION 

The  aim  of  this  investigation  is  to  contribute  to  the  understanding 
of  the  dynamics  of  large-scale  upper  oceanic  fronts,  with  particular 
emphasis  on  the  subarctic  front  of  the  North  Pacific.  Although  this 
study  is  primarily  intended  to  be  a  numerical  one,  theoretical  aspects 
such  as  the  stability  and  the  determination  of  the  time  scale  of 
dissipation  due  to  horizontal  friction  of  oceanic  fronts  also  have  been 
addressed. 

The  literature  related  to  oceanic  fronts  is  modest  while  the 
literature  of  atmospheric  fronts  is  quite  vast  and  profuse.  The  first 
attempts  in  the  understanding  and  the  physical  interpretation  of 
atmospheric  fronts  began  early  in  this  century.  The  classical  Norwegian 
school  was  the  first  which  attempted  to  understand  this  meteorological 
phenomenon.  The  papers  of  Bjerknes  (1919)  and  Bergeron  (1928)  are 
classical  in  this  matter. 

More  comprehensive  studies  of  atmospheric  fronts  showed  the 
inadequacy  of  the  classical  Norwegian  frontal  point  of  view.  The  study 
of  Sanders  (1955)  showed  the  evidence  of  a  new  type  of  front:  the 
surface  front.  The  intensity  of  this  frontal  zone  is  much  greater  at  the 
surface  than  at  higher  levels  in  terms  of  horizontal  gradients  of  both 
wind  and  temperature. 

The  papers  of  Reed  and  Sanders  (1953)  and  Reed  (1955)  draw  attention 
to  another  type  of  frontogenetlcal  phenomenon:  the  upper  tropospheric 
front. 
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The  papers  of  Hoskins  (1971),  Hoskins  (1972)  and  Hoskins  and 
3retherton  (1972)  should  be  interpreted  as  a  way  to  simulate  numerically 
these  two  fronts.  They  achieved  moderate  success  with  their  numerical 
simulations  by  making  numerous  assumptions:  Boussinesq  approximation; 
geostrophic  balance  across  but  not  along  the  front;  rigid  lid 
approximation;  omission  of  the  curvature  terms;  f-plane  approximation; 
and  unrealistic  deformation  field  which  cannot  be  valid  in  the  whole 
troposphere.  In  their  simulation  of  surface  fronts,  there  were  three 
salient  features:  intense  convergence  at  the  surface  front  and  a  strong 
upward  motion  above  it;  in  the  warmer  air  there  was  a  tendency  to  move 
upward  along  the  isentrops,  while  in  the  cold  air  there  was  a  tendency  to 
subside;  after  a  certain  time  the  numerical  solutions  showed  infinite 
vertical  vorticity  at  the  surface.  On  the  other  side  of  the  spectrum, 
their  simulation  of  the  upper  tropospheric  front  showed  a  remarkable 
discontinuity  between  tropospheric  and  stratospheric  air.  Hoskins  and 
Bretherton  were  successful  enough  in  simulating  the  folding  of  the 
tropopause,  but  only  down  to  400  mb.  Because  of  the  Boussinesq 
approximation,  the  vertical  displacements  were  limited.  Nevertheless, 
when  non-3oussinesq  effects  were  included,  Hoskins  obtained  only  slightly 
better  results.  Unfortunately ,  the  main  discrepancies  remain  unaltered. 

In  a  numerical  study,  Shapiro  (1975)  implemented  a  20-level 
isentropic  coordinate  primitive  equation  in  order  to  simulate  upper  level 
frontogenesis.  He  was  fairly  successful  in  predicting  the  folding  of  the 
trooopause  down  to  a  level  of  500  mb.  The  vertical  coordinate  he  used 
was  potential  temperature,  which  allows  a  vertical  variable  resolution. 

He  based  his  study  on  the  papers  of  Eliassen  and  Raustein  (1970),  using 
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their  same  initial  conditions.  In  our  model,  although  density  will  be 
the  vertical  coordinate,  the  philosophical  approach  will  be  somewhat 
similar  to  the  isentropic  coordinates  system  proposed  by  Eliassen  and 
Raustein  (1970).  Similar  studies  have  been  done  involving  research  with 
models  of  the  ocean,  rather  than  the  atmosphere. 

Tareyev  (1968),  using  a  simple  two-layered  model,  studied  the 
baroclinic  instability  in  the  section  of  the  Gulf  Stream  current  between 
Florida  and  Cape  Hatteras,  where  the  width  of  the  Gulf  Stream  current  is 

very  narrow  (approximately  100  km).  By  neglecting  the  3 -effect  and 

taking  into  accourt  the  non-geostrophic  effects,  he  was  able  to  show  that 

the  flow  is  barocl inical ly  stable  in  the  section  of  the  Gulf  Stream 

current  as  the  internal  Froude  number  is  less  than  one.  In  this  same 
paper  Tareyev  explained  the  contradiction  between  Kotchin's  model  and  the 
quasigeostrophic  models.  In  Kotchin's  model  the  baroclinic  instability 
occurs  when  the  vertical  velocity  shear  decreases,  while  in  the  latter 
models  the  baroclinic  instability  increases  with  increasing  vertical 
shear  after  reaching  a  critical  value. 

Orlanski  (1968)  applied  Kotchin's  model  in  an  f-plane  in  order  to 
study  the  instability  of  a  frontal  surface.  He  found  that  in  the  region 
of  Rossby  and  Richardson  numbers  less  than  three,  unstable  waves  exist 
for  all  wavelengths.  While  studying  the  region  of  complex  eigenvalues, 
he  found  a  new  form  of  instability  at  small  Richardson  numbers:  the 

♦ 

Rayleigh  instability.  Abramov,  et  al.  (1972)  extended  Orlanski 's  results 
for  a  3-plane  approximation.  Their  study  reveals  that  for  a  fixed  Rossby 
number,  baroclinic  instability  increases  rapidly  with  the  Richardson 

* 
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number  and  then  decreases  monotonically.  Moreover,  for  an  Increasing 
Richardson  number,  the  region  of  unstable  waves  narrows.  This  study  also 
reveals  the  destabilizing  effects  introduced  by  the  beta  effect  at 
certain  wavelengths.  In  all  quasi geostrophic  models  which  do  not 
consider  the  latitudinal  dependence  of  the  variables  ''f  the  basic  state, 
the  beta  effect  has  a  stabilizing  effect. 

The  study  and  understanding  of  oceanic  fronts  is  in  its  infancy  as 
compared  with  that  of  atmospheric  fronts.  There  are  two  powerful  reasons 
which  account  for  this  phenomenon.  While  an  enormous  amount  of 
meteorological  data  is  obtained  from  the  atmosphere  on  a  daily  basis,  it 
is  much  more  costly  and  difficult  to  obtain  accurate  data  from  the  ocean. 
Another  reason  is  that  in  the  atmosphere  at  a  certain  pressure  level,  the 
density  and  the  temperature  fields  are  inversely  related,  via  the 
equation  of  state. 

The  picture  is  completely  different  in  the  ocean.  Another  parameter 
is  present  in  the  equation  of  state:  salinity.  In  oceanography,  each 
field  will  have  to  be  studied  separately,  and  it  is  proper  to 
individualize  each  one  of  them.  Finally,  semantic  difficulties  arise 
from  the  fact  that  oceanic  fronts  exist  on  numerous  space  and  time 
scales.  Thus,  they  are  governed  by  different  dynamics,  characterized  by 
different  properties  and  driven  by  different  forces. 

When  the  temperature  and  salinity  fronts  happen  to  compensate  each 
other,  the  density  front  is  very  weak.  The  oceanic  subarctic  front 
typifies  this  case  (Roden,  1972).  In  an  opposite  sense,  in  the  case  of 
the  subtropical  front,  the  density  front  is  more  pronounced  because  the 
temperature  and  the  salinity  front  do  not  tend  to  compensate  each  other. 
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Hoskins  and  Bretherton  (1972)  noted  that  there  are  at  least  eight 
mechanisms  important  in  changing  temperature  gradients  and  forming 
atmospheric  fronts:  1)  horizontal  deformation  field,  2)  horizontal  shear 
mechanism,  3)  vertical  deformation  field,  4)  differential  vertical 
motion,  5)  latent  heat  release,  6)  surface  friction,  7)  mixing  processes, 
and  8)  radiation.  In  the  ocean,  the  first  six  mechanisms  can  be 
important;  the  last  two  are  negligible.  For  the  purpose  of  oceanic 
frontogenesis,  horizontally  non-uniform  buoyancy  fluxes  (heating  and 
cooling,  precipitation  and  evaporation,  river  run-off,  ice  melt,  ice 
brine,  etc.)  may  also  be  considered. 

31eck  (1978)  simulated  coastal  upwelling  frontogenesis  on  a  B -plane 
using  a  Lagrangian  type  of  formulation.  In  his  study  he  handled  the 
serious  problem  of  the  intersection  of  the  interface  with  the  sea  surface 
by  'using  a  new  technique:  the  "hole-filling  technique."  This  technique 
constitutes  no  more  nor  less  than  a  different  way  of  smoothing.  Clearly, 
"hole-filling"  is  not  the  answer  to  the  problem. 

Garvine  (1974)  developed  a  non-rotating  model  of  small-scale  oceanic 
fronts.  His  study  revealed  that  interfacial  friction  must  balance  the 
horizontal  pressure  produced  by  the  sloping  of  the  free  surface.  An 
intense  sinking  process  at  the  frontal  surface,  with  a  consequent  strong 
convergence  at  the  front  near  the  sea  surface,  was  oosarved.  He  extended 
this  study  and  developed  a  hyorodynamical  model  of  upper  ocean  fronts 
(Garvine,  1979).  Because  Garvine's  model  is  quasi-steady,  it  shows  the 
wind  stress  playing  a  secondary  role.  However,  it  is  well  known  that 
this  is  not  actually  the  case. 
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McVean  and  Woods  (1980)  considered  the  formation  of  a 
two-dimensional  upper  oceanic  front  forced  by  a  barotropic  horizontal 
stretching  deformation  field,  where  the  long  front  flow  is  geostrophic. 
While  assuming  a  deformation  rate  of  10"5  sec-1,  their  solutions  showed  a 
surface  density  discontinuity  in  three  to  four  days.  In  their  study,  the 
frontogenesis  in  the  temperature  field  depended  on  the  initial 
inclination  of  the  isotherms  to  the  isopycnics.  Furthermore,  the 
isotherm  pattern  penetrated  deeper  than  the  isopycnal  pattern. 

The  disparity  in  the  solutions  obtained  for  the  temperature  and  the 
density  fields  are  very  consistent  with  observations.  Typical  examples 
are  the  subarctic  front  and  the  subtropical  front,  in  the  winter  season, 
where  a  lack  of  barocl inicity  is  observed.  However,  because  of  the 
weakness  of  the  temperature  gradient  in  the  summer  season,  a  moderate 
baroclinicity  is  observed  in  the  subtropical  front  (Roden,  1974).  Cn  the 
other  hand,  both  the  Kuroshio  and  the  Gulf  Stream  fronts  are  strongly 
baroclinic  (Roden,  1976). 

Kao  (1980)  studied  the  dynamics  of  the  establishment  and  maintenance 
of  the  upper  oceanic  fronts  by  a  western  boundary  current.  The  dynamics 
are  controlled  by  two  dimensionless  parameters:  the  thermal  Rossby 
number,  Ro,  and  the  Ekman  number,  E.  In  his  study,  while  the  solutions 
proved  to  be  similar  for  a  large  Ro,  the  normalized  equations  depended 
strongly  on  E  alone.  After  a  quasi-steady  state  is  reached  and  for  E 
much  less  than  unity,  a  front  is  formed.  An  intense  geostrophic  jet 
along  the  front,  a  cross-front  circulation  with  surface  discharge  towards 
the  front,  and  a  return  flow  at  a  greater  depth  were  observed.  Potential 
vorticity  was  conserved  in  the  less  dense  waters. 
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In  this  study  a  numerical  model  has  been  developed  to  better 
understand  the  dynamics  of  the  subarctic  front.  Instead  of  following  the 
traditional  path  of  spinning  up  the  ocean  model  from  rest  in  order  to 
"build  up"  a  front,  a  geostrophically  balanced  oceanic  front  is  imposed 
as  an  initial  condition.  The  fronts  to  be  simulated  will  have  two 
different  forms.  In  the  first  case,  the  front  will  be  parallel  to  the 
latitudinal  circles,  while  in  the  second  instance  the  front  will  have  an 
undulatory  initial  form  around  42°N.  A  set  of  different  atmospheric  wind 
forcing,  representing  a  variety  of  different  atmospheric  conditions,  is 
imposed  upon  these  two  oceanic  fronts.  The  dynamic  response  of  these  two 
fronts  to  the  different  atmospheric  forcing  is  carefully  studied. 

Because  a  four-layer  numerical  model  is  employed,  a  better 
understanding  of  the  dynamical  response  of  the  ocean  at  the  wake  of  a 
storm  is  achieved.  As  a  matter  of  fact,  the  previous  studies  that  dealt 
with  the  cyclone's  problem  employed  a  two-layer  model.  Basically,  our 
results  show  a  close  agreement  with  Geisler’s  linear  solutions  (1970). 
This  research  confirms  that  water  is  upwelled  in  the  first  300  meters  at 
the  wake  of  the  cyclone.  This  upwelling  is  of  the  order  of  20  meters. 

Vertical  oscillations,  related  with  inertia  gravity  waves  induced  by 
the  storm's  passage  are  observed  (Chang  and  Anthes,  1979).  The  e-foiding 
time  scale  of  these  vertical  oscillations  is  calculated  to  be  of  the 
order  of  10  days.  This  value  confirms  the  fact  that  in  the  absence  of 
any  external  forcing,  the  induced  kinetic  energy  is  quickly 
redistributed.  After  a  week  of  the  storm's  passage,  the  water  remains 
upwelled  in  the  first  300  meters.  It  can  be  concluded  that  upwelling 
will  continue  for  several  weeks. 


A  region  of  relative  maximum  convergence  is  observed  in  the  lower 
layers  below  a  region  of  relative  maximum  upwelling  in  the  upper  layer, 


which  may  be  explained  by  continuity. 

Roden  demonstrated  that  frontogenesis  is  strongly  favored  when 
3xx/3y  is  strongly  positive.  Such  is  the  case  for  the  subarctic  front 
which  happens  to  be  located  to  the  south  of  the  maximum  westerlies, 
which,  in  turn,  will  generate  a  convergence  zone  in  that  region.  This 
particular  result  has  been  confirmed  during  the  course  of  the  present 
research  when  the  classical  case  of  westerlies  to  the  north  and 
easterlies  to  the  south  of  the  front  is  simulated. 

Parallel  to  the  numerical  aspect  of  the  present  research,  some 
theoretical  problems  are  also  addressed.  One  of  these  problems  is  to 
study  the  barotropic  stability  of  the  subarctic  front.  It  has  been 
proved  that  if  |~u-ci2>g-TT  TT'2  ap/Af^  everywhere,  the  front  is  stable. 
Another  theoretical  result  of  the  present  study  shows  that  the  time  scale 
of  dissipation  of  horizontal  friction  is  extremely  large. 

It  has  been  observed  that  numerical  instability  is  caused  when  the 
constancy  in  the  density  jump  between  layers  is  maintained  during  the 
whole  time  integration.  This  instability  is  usually  overcome  by  removing 
such  a  constancy  (McNider  and  O'Srien,  1973).  A  detailed  analysis 
explaining  the  reason  for  this  numerical  instability  has  been  conducted. 
In  decomposing  the  vertical  modes  of  the  linearized  model  equation,  it  is 
shown  that  the  phase  speed  of  the  baroclinic  modes  is  larger  in  the  case 
where  the  density  jump  is  constant  than  in  the  case  of  the  density  values 
arbitrarily  chosen. 


i 


9 


This  research  develops  a  better  understanding  of  frontal  oceanic 
response  to  atmospheric  variations.  With  a  simple  numerical  model,  the 
principal  atmospheric  forcings  due  to  the  wind  action  have  been 
simulated.  In  some  cases,  these  results  can  be  extended  easily  to  other 
oceanic  fronts  located  in  other  geographical  areas.  However,  there  are 
still  many  problems  that  need  to  be  studied.  Thermodynamic  processes 
such  as  radiative  fluxes,  precipitation  and  evaporation  need  to  be 
implemented  in  a  numerical  model  of  an  oceanic  front.  The  numerical 
simulation  of  such  phenomena  will  definitively  represent  a  new  challenge 
to  the  numerical  modeller.  A  different  track  will  have  to  be  followed 
without  the  present  advantage  of  using  density  as  a  vertical  coordinate, 
as  density  will  no  longer  be  constant  either  in  space  or  in  time.  In 
other  words,  a  hydrothermcdynamical  model  will  have  to  be  implemented  in 
such  a  way  that  the  simple  Lagrangian  formulation  of  the  numerical  model 
used  in  this  study  will  no  longer  be  valid. 


2.  THE  MODEL 

Consider  a  stably  stratified,  hydrostatic,  rotating  fluid  where  the 
Boussinesq  approximation  is  made.  The  model  consists  of  five 
incompressible  layers  which  have  initial  densities  pj  and  thickness  hj, 
so  that  j  =  1(1)5,  counting  down  from  the  upper  layer.  The  barotropic 
wave  modes  are  characterized  by  a  very  large  propagational  speed  compared 
with  that  of  the  baroclinic  wave  modes.  Because  the  lowest  layer,  layer 
5,  is  assumed  to  be  infinitely  deep,  we  can  consider  it  to  be  at  rest 
everywhere.  With  this  assumption,  all  the  external  barotropic  modes  are 
filtered  out.  The  free  surface  of  the  uppermost  layer,  5,  remains  at  a 
height  H  above  the  ground.  It  is  obvious  that  the  relation 

5 

H  -  [hj  (1) 

i=l 

must  hold.  The  j-th  layer  has  an  initial  constant  thickness,  hj,  a 
constant  density  pj,  and  a  pressure  Pj  of  the  form 

j-1  5 

Pj  =  Pa  +  +  g  [  p fhj  +  gpj  [  l  hj-z]  (2) 

i=l  1-j 

at  any  point  (x,  y,  z)  inside  this  layer,  where  P3  is  the  atmospheric 
pressure. 

Because  the  lowest  layer  is  considered  at  rest  during  the  who1e  time 
integration,  the  following  relation  holds: 

5 

7[  I  p ihi  +  Px5j  =  0  (3) 

i*l 

where  v  represents  the  horizontal  gradient. 
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Eliminating  vh5  and  75  from  (1)  and  (3)  yields 

5 

7h  =  -V  I  hi 
i=2 


5 

V5  =  -7  l  (P1/p  )hi 
i=l 


(4) 

(5) 


Neglecting  horizontal  variations  of  the  atmospheric  pressure,  Pa, 
be  rewritten  as 


j-1  5 

7Pj  =  g7[p  5  +  l  pjh,  +  pj  l  hi] 
i=l  i=j 


(2)  can 

(5) 


The  substitution  of  (4)  and  (5)  in  this  last  expression,  after  some 
algebraic  manipulations,  yields 


5  j-1  5 

vPj  =  -gv[  l  Pi  hi  -  l  Pi  hi  -  Pj-  l  hi]  (7) 

i=2  i=2  i = j 


Therefore,  the  pressure  gradient,  7Pj,  will  have  an  expression  of  the 


form 


vPj  =  gv[  I  (Pj-Pi)hi] 
i=j+l 


(3) 


Equation  (3)  is  not  very  useful  as  it  stands  because  it  contains  terms  of 
the  lowest  layer  (which  remains  motionless  during  the  whole  time 
integration).  Therefore,  it  is  advisable  to  substitute  h. .  The 
following  relation  holds: 

4 

7h_  =  -7  7  hj<  (9) 


The  substitution  of  (3)  into  (9)  yields  a  final  expression  for  7P.-  of  the 

J 


form 


7P  i 


4 

9  L  {Pc-Ov'/Th-; 

i  =  l  3 


(10) 


where  v  =  max  ( i ,  j ) . 
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Taking  into  consideration  equation  (10)  the  basic  horizontal 

momentum  equations  for  the  j-th  layer  may  be  expressed  as 

duj  4  p  _p  3 h-j  Xx 

—  -  fv-i  =  -g  l  [U-^il  —  +  Av2U]-  +  - 

dt  J  aiti  Pj  J  3x  J  Pi  hi 


dv 

dt 


3  hi 


i  ^o-o  0!,i  tY 

-  +  fuj  =  -g  l  [-5 — iij  —  +  AV2v-j  +  - 

^  J  Y=i  pj  J  ay  J  Pihi 


(11) 

(12) 


d  3  3  3 

where  f  =  f0  +  3y  and  —  =  —  +  us  —  +  vs  —  .  It  is  worthwhile  to 

u  dt  3t  J  3x  J  3y 

notice  that  the  wind  stress  forcing  is  exerted  only  on  the  top  layer. 

The  continuity  equation  for  the  j-th  layer  will  have  the  form 


3hj  3  3 

—  +  —  (uihi)  +  —  (vjhj)  =  0 
3t  3x  J  -  3y  J  J 


(13) 


A  Cartesian  coordinate  system,  x  and  y,  and  the  corresponding  velocity 
components,  u  and  v,  (which  are  positive  eastward  and  northward, 
respectively) ,  is  used.  The  vertical  coordinate  is  density. 

This  class  of  model  tends  to  be  numerically  unstable  when  the 
interface(s)  reaches  the  sea  surface.  A  closer  look  at  the  Droblem 
reveals  that  this  numerical  instability  is  most  commonly  caused  by  the 
imposition  of  closed  walls.  While  working  with  closed  boundaries,  the 
accumulation  or  removal  of  water  mass  at  both  northern  and  southern 
boundaries,  caused  by  the  Ekman  drift,  tends  to  increase  in  an 
unrealistic  form.  As  extreme  atmospheric  forcing  functions  are  applied 
in  this  study,  this  type  of  numerical  instability  appears  very  quickly. 
As  this  is  not  an  acceptable  solution  to  the  problem,  an  alternative  is 
required.  On  the  other  hand,  from  an  objective  point  of  view,  it  is 
philosophically  inconsistent  to  impose  closed  walls  where  they  do  not 
exist. 
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The  approach  taken  here  is  quite  simple.  By  scaling  the  x-momentum 
equation,  the  two  leading  terms  are  selected.  These  are  the  Coriolis  and 
the  wind  stress  terms.  Now  the  meridional  velocity,  v,  can  be  determined 
at  the  northern  and  southern  boundaries  as  a  function  of  the  wind  stress. 
As  a  consequence,  the  Ekman  drift  at  both  boundaries  is  prescribed 
instead  of  closed  walls  being  imposed.  This  new  approach  constitutes  a 
clear  improvement,  as  the  numerical  instability  caused  by  the 
intersection  of  the  interface  with  the  sea  surface  at  the  boundaries  is 
avoided.  Furthermore,  with  the  present  approach  there  is  no  constraint 
either  to  terminate  the  time  integration  or  to  switch  to  a 
hydrothermodynamical  system. 

The  model  geometry  is  an  idealization  of  the  North  Pacific  Ocean. 

At  both  eastern  and  western  boundaries,  periodic  boundary  conditions  are 
implemented.  The  east-west  extension  of  the  basin  is  approximately 
13,000  km,  which  roughly  represents  the  distance  between  the  two 
continental  areas.  Of  particular  interest  is  the  dynamics  of  the 
mid-ocean,  in  particular  in  a  subregion  of  1,000  km  along  the  x-axis.  A 
stretching  coordinate  is  implemented  along  this  axis  in  such  a  way  that 
the  number  of  grid  points  close  to  the  eastern  and  western  boundaries  is 
minimal.  At  the  same  time,  the  distance  between  grid  points  is  constant 
in  the  region  of  interest. 

The  same  philosophy  is  applied  for  the  meridional  coordinate. 

Because  a  better  resolution  is  desired  in  the  middle  of  the  N-S  extent  of 
the  basin,  where  the  oceanic  front  is  initially  centered,  a  stretching 


coordinate  is  employed  along  the  y-axis  in  a  way  such  that  the  distance 
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between  grid  points  near  the  front  is  of  the  order  of  4  kms.  Since  there 
is  no  interest  in  considering  the  boundary  effects  at  the  northern  and 
southern  boundaries,  the  basin  is  extended  from  30-54°N,  while  the  region 
of  interest  lies  between  40  and  44°N.  Again,  the  number  of  grid  points 
away  from  the  region  of  interest  is  reduced  to  a  minimum.  Therefore,  the 
region  of  interest  consists  of  a  flat-bottomed  ocean  basin,  whose 
east-west  extension  is  on  the  order  of  1,000  kms,  while  the  north-south 
extension  is  on  the  order  of  450  kms. 

Because  of  the  particular  interest  in  considering  the  dynamical 
response  of  the  subarctic  front  caused  by  a  cyclone  travelling  from  north 
to  south,  a  second  rectangular  basin,  following  the  same  idea  described 
above,  is  implemented.  In  this  case,  the  inner  mesh  has  an  extension  of 
200  kms  and  1,000  kms  in  the  E-W  and  N-S  directions,  respecti vely.  On 
the  other  hand,  the  dimensions  of  the  outer  mesh  in  the  x-  and 
y-directions  are  of  1,000  kms  and  5,900  kms,  respecti vely. 

Two  initial  formulations  of  oceanic  fronts  are  considered.  In  both 
cases  the  oceanic  fronts,  centered  in  the  middle  of  the  north-south 
extent  of  the  basin,  are  determined  by  a  hyperbolic  tangent  of  the  form 

y*-W/2 

hj  =  hj0  -  Ahj  tanh( — - - )  (14) 

Ly 

where  the  subindex  j  represents  the  j*th  layer  so  that  j  =  1(1)5, 
counting  down  from  the  upper  layer  and  W  represents  the  width  of  the 
basin.  In  the  first  formulation 

y*  3  y»  (15. a) 

while  in  the  second  one 


y*  =  y  A  3  sin  kx  , 


(15. b) 
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where  B  (  *  Ly/3)  is  the  amplitude  and  k  (  =  2n /X )  is  the  zonal  wave- 
number.  The  value  of  the  wavelength,  X,  is  arbitrarily  chosen  to  be  500 
kilometers,  so  as  to  model  an  initially  meandering  front.  These  oceanic 
fronts  are  in  geostrophic  balance  with  the  Coriolis  parameters  in  such  a 
way  that 

4  n  -o  3h-j 

uj  3  -(g/f)  I  C-5— ^)r—  d6) 

i*i  pj  3y 

By  continuity,  the  meridional  velocity,  v,  in  the  first  case  (front 
without  waves)  is  equal  to  a  constant,  which  is  set  to  be  equal  to  zero. 
In  the  second  formulation,  because  of  the  undulatory  initial  form  of  the 
oceanic  front  around  42°N,  the  meridional  velocity  is  determined  by 

4  p  -o  3  hi 

Vj  -  (g/f)  I  (-5-*)—  (17) 

i  =  l  Pj  3  x 

It  is  desirable  to  study  the  dynamical  response  of  the  temperature 
and  the  salinity  fronts  at  the  wake  of  a  strong  cyclone.  For  this 
purpose  an  initial  temperature  front  is  implemented  in  the  model  via  an 
hyperbolic  tangent  of  the  form: 

y*- W/2 

Tj  =  T.j  -  ATj  tanh(  — - — )  (18) 

y 

where  Tj  represents  the  temperature  of  the  j*th  layer  (Fig.  la).  An 
initial  salinity  front  is  simulated,  via  the  equation  of  state,  by  the 
relation: 

Sj  =  (pj  -  Xj°  -  arj)/b  (19) 

where  Sj  represents  the  salinity  of  the  j*th  layer,  and  a  and  b 
represents  the  coefficients  of  thermal  expansion  and  saline  contraction, 
respectively  (Fig.  lb).  The  constant  Xj°  has  the  form: 
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Xj°  =  p°  +  a^3  -  bS°  (20) 

where  the  superscript  "°"  represents  a  state  of  reference.  The  model 
representation  of  the  temperature  and  salinity  fronts  is  in  good 
agreement  with  Roden's  (1976)  observations  of  such  fronts. 

The  temperature  is  prognostically  evaluated,  at  time  step  n  +  1,  by 
the  relation: 


d? 

dt 


0 


(21) 


while  the  salinity  is  diagnostically  evaluated  by  (19). 

No  flux  of  temperature  or  salinity  across  the  boundaries  is  allowed. 
IT  3  S 

Thus,  the  relation  —  =  —  =  0,  holds  for  the  four  boundaries,  where  the 
3  3n  3n 

symbol  —  stands  for  the  normal  derivative. 

3n 

The  atmospheric  forcing  used  in  this  research  is  rather  extreme. 
Although  the  oceanic  fronts  are  imposed  as  an  initial  condition,  the 
severity  of  the  different  atmospheric  forcing  being  used  may  be  justified 
by  the  fact  that  this  research  focuses  on  the  dynamical  effects  that 
these  different  forcings  exert  on  the  oceanic  front  in  a  very  short  time 
scale.  A  more  practical  reason  is  the  purpose  of  saving  computer  time 
during  the  numerous  try-outs. 

One  case  study  is  the  dynamical  effects  caused  by  the  atmospheric 
westerly  winds.  For  this  purpose,  a  wind  stress  of  10  dynes  cm"2  blows 
toward  the  eastern  boundary  at  all  grid  points  during  the  whole  time 
integration.  Another  case  study  is  the  effect  of  the  passage  of  a  severe 
atmospheric  front.  Initially  located  at  1,000  km  to  the  west  of  our 
region  of  interest,  the  atmospheric  front  travels  across  the  basin  at  a 
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speed  of  500  kms/day  from  the  NW  to  the  SE  corner.  Behind  the 
atmospheric  front,  the  northwesterly  winds  are  represented  by  the 
formulation 

xy  -  -T0  sine  (22. a) 

t*  *  T0  cose  (22.b) 

where  the  amplitude  of  the  wind  stress,  T0,  and  the  angle  that  the 
atmospheric  front  forms  with  the  x  axis,  9,  are  given  arbitrary  values  of 
7  dynes/ cm2  and  45°,  respecti vely. 

The  effects  caused  by  the  passage  of  a  strong  atmospheric  cyclone 
are  considered.  A  westerly  wind  stress  of  2  dynes/cm2,  initially 
imposed,  is  suddenly  interrupted  by  the  passage  of  this  cyclone  traveling 
towards  the  eastern  boundary  at  a  speed  of  500  km/day.  The 
circumferential  velocity  of  this  cyclone,  u,  is  formulated  in  the 
following  manner: 

Vo(l+r/ra)“l  10<rs500  kms 

« •  m  *  l  (23) 

0  r<10  kms 

where  r2  =  rx2  +  ry2  ,  and  v0  is  equal  to  70  m.  s"1 . 

Having  determined  u,  the  wind  stress  of  the  cyclone  may  be  computed 

as  follows 

tX  3  Pa  cD|V!  u  (24. a) 

'■'J  -  Pa  cD {V |  v  (24.b) 

Typical  values  of  pa  (  *  10“3  g  cm"3),  cq  (  =*  2.5.10"3)  and  r0 

(  =»  75  km)  are  used. 

The  classical  effect  of  westerlies  at  the  north  and  easterlies  at 
the  south  of  the  oceanic  front  is  also  addressed.  The  curl  is 
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represented  by  a  cosine  function 

tx  =  -T^  cos  (fry*/W)  (25) 

where  T  is  given  the  value  of  10  dynes  cm-2.  The  value  of  x 'J  is  set  to 
be  equal  to  zero. 

For  the  purpose  of  saving  computer  time  and  in  order  to  avoid  the 
numerical  noise  produced  by  the  higher  wave  number,  the  model  equations 
are  integrated  in  a  staggered  grid  in  space.  The  staggering  of  the  grid 
is  such  that  the  zonal  velocity  points,  u,  alternate  with  height  points, 
h,  along  the  latitudinal  direction,  while  the  meridional  velocity  points, 
v,  alternate  with  height  points,  h,  along  the  meridional  direction. 

The  nonlinear  terms  are  expressed  in  the  "semimomentum"  formulation. 
This  scheme  is  used  to  ensure  the  elimination  of  fictitious  sources  of 
kinetic  energy. 

The  model  consists  of  four  active  layers  of  depth 

^  =  3K,  h4  =  h3  =  h£  ■  K,  K  =  50  m.  (26) 

The  first  layer  is  chosen  to  be  that  deep  because  it  has  been  observed 
that  both  the  salinity  and  temperature  fronts  tend  to  compensate  in  the 
first  150  meters. 

The  vertical  density  profile  is  constructed  so  that  there  is  a 
difference  of  2.10-3  g  cm-3  between  the  first  and  second  layers,  a 
difference  of  3. 10 — *  g  cm"3  between  the  two  bottom  layers,  and 
differences  of  5. 10“^  g  cm-3  between  the  remaining  layers. 


3.  ON  THE  ROLE  OF  HORIZONTAL  FRICTION 


An  attempt  to  determine  the  time  scale  of  dissipation  of  an  oceanic 
front  is  made.  The  study  is  conducted  in  such  a  way  that  if  the  time 
scale  of  the  dissipation  of  the  oceanic  front  is  larger  than  the  time 
scale  of  the  perturbation  produced  by  the  wind  stress,  it  is  possible  to 
neglect  lateral  friction.  Because  the  longitudinal  frontal  length  scale, 
Ly,  is  smaller  than  the  latitudinal  frontal  length  scale,  Lx>  the 
derivatives  with  respect  to  x  will  be  neglected  as  compared  to  the 
derivatives  with  respect  to  y.  A  simple,  two-layer  model  is  considered. 
The  lower  layer,  layer  2,  is  chosen  to  be  infinitely  deep  in  such  a  way 
that  the  barotropic  mode  is  filtered  out.  The  trivial  solution 
u2  =  v2  =  0  follows  quite  naturally. 

The  pressure  in  the  upper  layer,  layer  1,  will  have  the  form 

Px  =  Pa  +  p^CV2)  +  p23(H”hi)  (27) 

where  Pa  represents  the  atmospheric  pressure  and  z  is  measured  from  the 
bottom.  The  gradient  of  PL  yields 

7Pl  =  -pl  g*7hx 


(28) 


P  -P 

where  g*  =  g (-1—2}  represents  the  reduced  gravity.  The  horizontal 
Pi 

fluctuations  of  the  atmospheric  pressure  are  neglected  in  the  present 
study  of  friction.  For  convenience,  the  subindex  1,  representing  the 
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upper  layer,  will  be  omitted.  The  equations  of  motion  for  the  upper 
1 ayer  are 


3u 

n 


+ 


3U 

V  — 

3y 


(29) 


3  V  3v  3h  32v 

—  +  v—  +  fu  =  -g*—  +  A— - 
3 1  3y  3y  Sy2 


(30) 


3h  3 
n  3y 


(vh)  =  0 


(31) 


where  the  forcing  due  to  the  wind  stress  terms  is  neglected  in  order  to 
retain  only  the  forcing  due  to  the  frictional  terms. 

This  non-linear  set  of  equations  is  not  easy  to  solve.  It  is 
assumed  that  the  frictional  effect  will  slowly  change  a  steady,  inviscid 
initial  front,  and  a  perturbation  solution  is  sought  by  assuming  that 
friction  has  a  weak  effect. 

For  this  steady  oceanic  front  to  persist,  the  interface,  h0,  must  be 
balanced  initially  by  a  zonal  velocity  of  the  form 

q*  *  hg 

u°  ■  rr  (32) 

where  the  subscript  "o'*  stands  for  the  basic  state. 

The  interface  has  a  mean  value,  h00,  while  its  variation  is  chosen 
to  have  a  hyperbolic  tangent  centered  in  the  middle  of  the  N* S  extant  of 
the  basin.  Therefore,  the  initial  thickness  of  the  upper  layer  is 
represented  as 

,y- w/2. 

■h0  ■  hg0  -  oh  tanh  i  — — j  (33) 

Ly 

where  oh  represents  the  maximum  amplitude  of  the  interface  perturbation. 
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8y  continuity,  the  meridional  velocity,  v,  is  initially  set  to  be  zero. 

The  viscous  perturbation  solutions  are  obtained  by  assuming  the 
value  of  A  to  be  very  small.  Let 


u  =  u0  +  u1  (34a) 
v  s  v0  +  v'  (34b) 
h  =  ho  +  h'  (34c) 


where  the  perturbation  variables  (which  are  of  order  A)  are  denoted  by 
primes.  The  equations  for  the  upper  layer  then  become 


3  U ' 

at" 


a2h0 

ay2 


-fv' 


,Ag*  33flo  32 u ' 


(35) 


3  v' 

at" 


3  V  1 

+  >/' -  +  f  u0  +  f  u 1 

ay 


ay2 


(36) 


ah' 

a"T 


a 

ay 


(V) 


(v'h1 ) 


0 


(37) 


a  u'  a  v 1  a2u‘ 

The  terms  v' — ,  v' — ,  a(v'h')/3y,  A — — ,  and  A32v'/ay2  are  of 
ay  ay  ay2 

second  order.  Tnus,  they  are  negligible.  The  terms  fu0  and  -g*ah0/3y 

represent  the  basic  state  balance.  In  order  to  filter  out  the  gravity 

waves,  the  term  av'/at  is  omitted.  With  these  simplifications,  the 

above  equations  reduce  to  linear  perturbation  equations  of  the  form 

3u' 

at” 


g*v'  a2h0 

— ♦  */«*)  * 

, Ag*,a3  ^o 

"f”J  a” 

(33) 

u'  =  -(g*/f)  — 

ay 

(39) 

ah'  a 

~  +  ay  (v  h°'  3  0 

UO) 
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To  solve  this  system,  a  single  equation  for  v‘  is  first  needed. 
Replacing  the  value  of  u'  in  (38)  yields 

a3  h0 


g*  32h  g*  a2h0  Ag* 

-  t _ L  .  L  v'f _ -  +  f2/g*l  =  -f-L] - 

f  3y3t  f  l3y2  73  J  '  f  sy3 

3 

Applying  the  operator  ( g*/f ) —  to  (40),  it  follows  that 

3y 

q*  a2  h  *  32 

~)( -  +  -7  (v'h0))  =  0 

f  3t3y  3y2  u 

The  addition  of  (41)  and  (42)  yields 

32  V  ,a2h0/3y2+f2/g*^  33hn 


(41) 


(42) 


ay2 


;V  •  -A77- 
3y3 


(43) 


where  V  (=  v'h0)  is  the  mass  transport  in  the  upper  layer.  The  boundary 
conditions  are  such  that  the  meridiona7  mass  transport,  V,  is  set  to  be 
zero  at  both  the  northern  an  southern  boundaries.  It  is  of  interest  to 
note  that  v  is  independent  of  time. 

As  h0(y)  is  a  known  function,  the  coefficients  of  (43)  are  known  and 
this  equation  can  be  rewritten  as 


32  V 

-r  -  F(y)V  »  -G(y) 
ay2 


where 


and 


r(y) 


32  h( 

3~ 


+  f2/g* 


s(y) 


ho 

a3h0 

A_ 

oy3 


(44) 


(45a) 


(45b) 


In  order  to  cover  a  wider, range  of  oceanic  fronts,  two  different 
meridional  length  scales,  Ly,  are  considered.  The  chosen  dimensions  are 
1  and  50  kilometers,  respecti vely.  In  the  first  case,  an  approximate 


l 
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value  of  F(y)  is: 


F(y) 


.  fVg* 


(46) 


while  in  the  second  case,  F(y)  has  an  approximate  value  of: 

.  32h0/3y2 

F(y) - 


(47) 


Using  the  first  approximation,  (44)  may  be  rewritten  as: 

32V  33h0 

.  -  a2V  =  -A - 

3y2  3y3 


(48) 


where  a2  =  f2(g*h0)-1  is  the  inverse  of  the  Rossby  radius  of  deformation 
squared.  For  this  particular  case,  an  analytical  solution  can  be 
obtained.  If  the  Ekman  number  is  set  to  be  equal  to  0.1,  the  resulting 
value  of  A  is  then  10  rr^sec"1.  The  scaling  of  (48)  shows  that  the  first 
and  third  terms  are  the  largest  ones,  by  two  orders  of  magnitude. 
Therefore,  (43)  has  the  form 


32V 

a"? 


2AAh 


/  , ,y-W/2,  \ 

1 2sinh2(-— - — j  -  1\ 

y 


\cosh4(~^) 


The  analytical  solution  of  this  expression  is 


(49) 


,,,  .  -^h  y- W/2, 

/(y)  =  -  sac.r-' - 

Ly  Ly  ' 


'50) 


7ne  integration,  with  respect  to  time,  of  the  continuity  equation 


yields  a  perturbation  expression  for  the  interface,  h',  of  the  form 


3V 

h  3  —  t  (51) 

3y 


U 
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Considering  the  worst  possible  situation,  the  smallest  time  scale  of 
dissipation,  T,  is  defined  in  such  a  way  that 

3  V  "l 

T-h'|jjlM  (52) 

3 

where  I~!m  represents  the  maximum  absolute  value  of  such  derivative. 

The  value  of  h'  is  arbitrarily  set  up  to  be  equal  to  10  meters.  For 

,3V, 

different  values  of  A  and  L,  different  values  of  |— im  may  he  obtained. 

A  list  summarizing  these  cases  is  condensed  in  table  1. 

However,  this  solution  is  restricted  to  oceanic  fronts  of  the  order 
of  ore  kilometer.  As  it  is  stated  previously,  the  subarctic  front 
has  a  meridional  length  scale  of  the  order  of  thirty  kilometers.  A 
new  scaling  of  (43)  using  this  length  scale,  reveals  that  the  three  terms 
of  this  equation  are  comparable.  Thus,  the  full  equation  has  to  be 
solved. 

In  using  a  meridional  length  scale  of  50  kilometers,  (48)  has  the 

form: 

3  2V  s2ho  V  33ho 

3  y2  3^  h0  A5y3 

The  general  solution  of  the  corresponding  homogeneous  equation  (53), 
Vc(y),  is: 

Vc(y)  =  C ]_ V !  -  C2V2  (54) 

where  'J\  =  h0,  V2  =  'n0J’h0"2dy  and  Cj  and  C?  are  arbitrary  constants.  The 
method  of  variation  of  parameters  is  used  to  solve  the  nonhonogeneous 
equation  (48).  This  method  requires  the  replacement  of  the  constants  C]_ 
and  C2  by  two  arbitrary  functions,  x(y)  and  c(y).  Tnese  arbitrary 
functions  will  be  determined  in  such  a  way  that  the  particular  solution. 


27 


Vp(y)»  has  the  form 

vp(y)  -  x(y)Vi(y)  +  t(y)  v2(y)  (55) 

The  arbitrary  function,  x  and  c,  must  satisfy  certain  conditions  in 
order  to  satisfy  the  nonhomogeneous  O.D.E.  These  conditions  are 


Xyvl  +  Syv2  =  0 

XyVl y  +  ?yV2y  =  V(y)  (56) 

where  xy»  ?y»  Viy  and  V2y  represent  the  meridional  derivatives  of  the 
functions  x»  Vj  and  V2  respectively ,  and  v (y)  represents  the  forcing 


function  defined  by 


f(y)  =  r 


2AhA^ 

y 


2s i nh2  C  — “3  -  l" 


y-W/2, 
cosir  1 - j 

'  Ly 


(57) 


The  system  of  equation  (56)  gives  us  a  solution  for  the  functions 
Xy  and  Cy.  These  values  are: 

Xy  «  -V2(y)  *(y)/W(Vi,  V2)  (53a) 

Cy  =  Vi(y)  ?(y)/W(V1,  V2)  (58b) 

where  the  wronskian,  W(V;,  V2),  is  equal  to  one. 

The  values  of  x  and  z,  at  each  grid  point,  are  determined  by 
numerical  integration.  To  achieve  this  purpose,  Simpson's  rule  is  used. 
Knowing  the  values  of  x  and  z,  a  final  expression  for  the  meridional  mass 
transport,  V,  can  be  numerically  evaluated  at  each  grid  point.  Such  an 
expression  has  the  form 


v(y)  *  (Ci  +  x)  vi(y)  +  (C2  +  O  v2(y)  (59) 


dV(y) 

Again,  by  plotting  the  function  -  versus  y,  a  maximum  absolute  value 

dy 

of  the  former  function  may  be  determined  for  different  values  of  A  and 


Ly.  Table  2  summarizes  these  results. 
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The  analysis  of  Tables  1  and  2  shows  the  same  salient  features: 
a)  The  larger  the  coefficient  of  diffusion.  A,  keeping  Ly  constant,  the 
smaller  the  time  scale  of  dissipation,  T;  b)  The  larger  the  meridional 
length  scale  of  the  oceanic  front,  Ly,  keeping  A  constant,  the  larger  the 
time  scale  of  dissipation. 

From  observations,  it  is  well  known  that  the  meridional  length  scale 
of  the  subarctic  front  is  of  the  order  of  thirty  kilometers.  Thus,  the 
results  presented  in  Table  2  take  full  account  of  the  different 
alternatives  that  might  occur.  These  results  reflect  the  same  basic 
conclusion:  friction  plays  no  significant  role,  as  the  time  scale  of 
dissipation  due  to  horizontal  friction  is  very  large. 
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Table  2. 

Same  as  Table  1, 

but  using  equation 

(53). 

4.  ANALYSIS  OF  THE  NUMERICAL  EXPERIMENTS 

There  is  a  considerable  lack  of  numerical  simulations  in  the 
modelling  of  the  mid-oceanic  fronts.  Even  though  the  observational  and 
theoretical  studies  performed  by  Roden  are  useful  because  they  give  us  a 
fairly  good  idea  of  the  dynamics  of  the  oceanic  fronts  in  the  North 
Pacific  Ocean,  the  supporting  numerical  modelling  of  such  fronts  is 
1 acking. 

The  aim  of  this  research  is  to  fill,  in  part,  such  a  vacuum.  The 
emphasis  of  this  study  is  to  understand  the  role  the  atmospheric  forcing 
plays  on  the  dynamical  structure  of  the  subarctic  front.  A  complete 
series  of  different  atmospheric  wind  forcings,  which  account  for  all 
possible  atmospheric  situations,  is  designed  in  order  to  achieve  such  a 
goal . 

The  approach  used  here  is  similar  to  the  one  taken  by  Eliassen  and 
Raustein  (1970),  the  difference  being  that,  while  studying  the  atmos¬ 
pheric  front,  they  chose  potential  temperature  as  a  vertical  coordinate. 

In  this  research,  two  fronts  are  simulated.  Both  of  them  have  a 
different  initial  configuration.  The  philosophy  behind  the  second  choice 
is  that  the  meandering  of  oceanic  fronts  is  not  an  infrequent  phenomenon 
(Fig.  2). 

The  effects  of  the  atmospheric  wind  forcing  have  a  time  scale  of  a 
week  or  two  (Roden,  1977).  It  should  be  no  surprise  that  we  carry  out 
our  time  integration  in  sucn  a  time  scale. 
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KILOMETERS 


Figure  2.  Initial  configuration  of  the  height  field  of  the  subarctic 
front  when  y*  3  y  +■  3s  ink*  for  the  first  layer  in  the  (x,  y) 
plane.  The  values  in  the  figures  are  expressed  in  the  c.g.s. 
system. 
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Two  main  reasons  determine  our  choice  of  150  meters  for  the  depth 
of  the  first  layer.  On  one  hand  it  is  observed  that  within  that  layer 
both  the  temperature  and  the  salinity  fronts  do  compensate  each  other. 

A  very  weak  density  front  is  thus  observed.  On  the  other  hand,  a  high 
stability  layer,  whicn  may  be  used  as  an  index  of  the  depth  of  the  upper 
layer,  is  found  to  be  between  100  and  150  meters.  The  atmospheric 
forcing  is  limited  to  the  layer  between  the  sea  surface  and  the  high 
stability  layer  (Roden,  1970).  By  contrast  in  the  layer  between  150  and 
300  meters  the  density  front  is  pronounced,  as  the  temperature  and 
salinity  fronts  do  not  happen  to  compensate  each  other.  Our  experiments 
tend  to  confirm  the  fact  that,  in  the  case  of  the  subarctic  front,  the 
response  to  the  atmospheric  wind  forcing  is  confined  to  the  first  150 
meters.  The  simulation  of  a  strong  cyclone  constitutes  the  only 
exception. 

It  is  widely  known  that  the  subarctic  front  intensifies  or  decays  in 
response  to  atmospheric  forcings  (Roden,  1576).  Four  different  types  of 
atmospheric  forcings  can  be  identified:  momentum,  radiative  heat  fluxes, 
precipitation  and  evaporation,  and  turbulent  energy  fluxes.  Because 
frontogenesis  in  the  Central  Pacific  is  largely  dependent  upon  the 
configuration  of  the  wind  stress  field,  this  research  focuses  on  this 
particular  type  of  atmospheric  forcing. 

As  our  study  focuses  on  the  dynamics  of  the  mid  ocean,  the  lateral 
boundary  effects  are  not  our  concern.  The  task,  therefore,  is  to  set  up 
a  model  in  such  a  way  that  the  computational  boundary  effects  should  not 
affect  nor  contaminate  the  interior  solution.  To  achieve  this  thrust, 
the  model  is  constructed  so  that  the  boundaries  are  placed  far  apart 
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from  the  region  of  interest.  While  pursuing  this  goal,  we  try  to 
minimize  the  number  of  grid  points  outside  the  region  of  interest  by 
using  stretching  coordinates.  Our  first  analysis  is  over  an  area  of  1000 
kms  in  the  east-west  direction  and  450  kms  in  the  north-south  direction. 

The  approach  taken  here  is  that  a  geostrophically  balanced  oceanic 
front  is  imposed  as  an  initial  condition,  instead  of  the  traditional  way 
of  spinning  up  the  ocean  model  from  rest  in  order  to  "build  up”  an 
oceanic  front.  Our  approach  is  more  efficient  and  more  suitable  for  the 
purpose  of  this  research  as  the  ocean  model  can  be  run  for  a  short  time 
scale.  Obviously,  a  great  amount  of  computer  time  is  saved  using  this 
type  of  initialization. 

All  the  numerical  simulations  dealing  with  the  ocean's  response  to 
the  passage  of  a  strong  cyclone  had  been  performed  using  two  layer 
models.  A  deeper  understanding  of  the  physics  of  the  ocean's  response 
will  be  gained  with  the  implementation  of  multilayer  models.  This  is  the 
main  reason  why  a  four-layer  model  is  implemented  in  this  study.  For  a 
clearer  visualization  of  the  physics  irvulved,  the  pictures  (of  the  u,  v 
and  h  field  solutions)  presented  hereafter  will  reproduce  the  isoanomal 
field  only.  That  is  to  say,  the  departure  from  the  initial  state  will  be 
plotted. 

The  first  case  to  be  considered  is  the  most  frequent  one.  A  uniform 
westerly  wind  is  imposed  at  all  grid  points  during  the  whole  period  of 
integration,  10  days.  This  particular  time  scale  is  sufficiently  large 
to  give  a  good  idea  of  the  resulting  dynamics  of  the  oceanic  front  in 
response  to  this  wind.  The  Ekman  flow  to  the  right  of  the  westerlies 
causes  a  southward  displacement  of  the  subarctic  front  (Fig.  3).  At  the 
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Figure  3.  Isoanomalies  of  the  height  field  for  the  westerly  wind  case 
after  10  days.  The  values  in  the  figures  are  expressed  in 
the  c.g.s.  system. 
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same  time  a  region  of  upwelling  is  detected  at  the  site  of  the  front. 

This  upwelling  of  water  may  be  explained  by  the  southward  Ekman  flow. 
While  the  geostrophic  zonal  velocity  field,  u,  shows  no  change  at  all, 
the  meridional  velocity  field,  v,  shows  a  southward  component,  which 
accounts  for  the  southward  Ekman  transport  and  subsequent  downwelling  at 
that  boundary.  No  meanders  are  observed  along  the  front.  The  stability 
of  the  front  is  thus  verified. 

It  is  observed  that  frontogenesis  is  stronger  where  the  curl  of  the 
wind  stress  is  strongly  negative.  This  feature  is  observed,  in  the  case 
of  subarctic  frontogenesis,  on  the  southern  side  of  the  westerly  jet.  A 
convergence  zone  is  originated  as  a  consequence  of  this  phenomenon. 

The  classical  case  of  westerlies  at  the  north  and  easterlies  at  the 
south  of  the  front  is  also  considered.  Convergence  all  over  the  basin 
can  be  detected.  The  downwelling  observed  supports  the  observations  of 
Roden  (1977)  as  the  convergence  observed  favors  the  frontogeneticai 
effects  (Fig.  4). 

The  response  of  a  horizontally  homogeneous  ocean  to  the  passage  of  a 
strong  cyclone  has  been  studied  from  the  observational,  theoretical  and 
numerical  points  of  view.  Perhaps  the  most  important  observational 
research  in  this  context  was  performed  by  Leipper  (1967).  His 
observations  of  the  hurricane  Hilda  of  1964  in  the  Gulf  of  Mexico  showed 
the  evidence  of  water  having  been  upwelled  from  a  depth  of  60  m.  This 
kind  of  upwelling  explained  the  drop  of  5°C  in  the  sea  surface 
temperature  that  Leipper  recorded  in  the  lee  of  hurricane  Hilda.  7ne 
numerical  simulation  of  O'Brien  and  Reid  (1967),  tends  to  support  this 
evidence.  Implementing  a  two  layer  nonlinear  model,  these  authors 
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studied  the  ocean's  response  to  a  stationary,  axisymmetric  hurricane. 

They  were  able  to  detect  a  maximum  upwelling  of  30  meters  at  30 
kilometers  of  the  center  of  the  hurricane.  The  fact  that  upwelling  is 
responsible  for  the  lowering  of  the  sea  surface  temperature  is  in  sharp 
contradiction  with  the  basic  assumption  of  Jordan  (1964).  He  estimated 
that  mixing  processes,  such  as  precipitation,  evaporation  and  sensible 
and  turbulent  heat  transfer  were  responsible  for  such  a  phenomenon. 

The  observational  study  of  Black  and  Withee  (1976)  showed  the 
evidence  of  inertia-gravity  waves  in  the  wake  of  hurricane  Eloise  of 
1975.  Because  O'Brien  and  Reid's  (1967)  hurricane  model  is  stationary, 
it  could  not  account  for  these  storm  induced  inertia-gravity  waves.  In  a 
theoretical  study,  Geisler  (1970)  studied  the  linear  response  of  a 
two-layer  ocean  model  to  a  moving  hurricane.  Whenever  the  translation 
speed  of  the  cyclone  exceeded  the  baroclinic  long-wave  speed,  he  was  able 
to  detect  an  internal  wake  in  the  lee  side  of  the  storm. 

Chang  and  Anthes  (1978)  employed  a  numerical  model  to  study  the 
ocean's  response  to  a  moving  hurricane.  They  showed  that  vertical 
motions  associated  with  the  inertia-gravity  waves  persist  at  the  wake  of 
the  cyclone's  path.  Although  their  numerical  model  was  basically 
non-linear,  their  results  agree  with  Geisler's  linear  solutions. 

The  next  experiment  to  be  simulated  is  the  passage  of  a  strong 
cyclone  moving  from  west  to  east,  at  a  speed,  S,  of  50G  kilometers  per 
day,  along  the  centerline  of  the  front.  Initially  located  at  2500  km  to 
the  west  of  our  region  of  interest,  the  center  of  the  cyclone  reaches  the 
western  boundary  at  the  end  of  the  fifth  day.  At  the  end  cf  the  sixth 
day,  the  center  of  the  cyclone  is  located  in  the  middle  of  the  tone  of 
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interest.  Figures  5a  and  b  show  a  shift  to  the  south  of  the  maximum 

j 

current,  as  observed  by  Chang  and  Anthes  (1978).  ] 

The  center  of  the  cyclone  at  the  end  of  the  seventh  day  is  located  ! 

at  the  eastern  boundary  of  our  region  of  interest.  Figure  5c  shows  that 
the  region  of  maximum  uplifting  lags  the  center  of  the  storm  by  150 
kilometers  and  is  slightly  to  the  right  of  the  center  of  the  cyclone's 
path.  A  second  maximum  uplifting  is  located  about  500  kilometers  behind 
the  storm.  This  second  maximum  is  clearly  associated  with 
inertia-gravity  waves  at  the  wake  of  the  cyclone.  The  wave  that  is 
observed  in  the  lee  of  the  storm  has  a  wavelength  of  approximately  350 
kilometers.  This  value  is  very  close  to  the  one  observed  by  Black 
(1977),  450  kilometers,  along  the  wake  of  hurricane  Ginger  of  1971,  and 
the  one  obtained  by  Chang  and  Anthes(1979) ,  of  450  kilometers.  The  main 
differences  between  our  model  and  Geisler's  is  that  our  model  is 
nonlinear  and  that  the  translation  of  the  cyclone  is  slower  in  our 
simulation.  These  two  factors  may  account  for  the  slight  discrepancy  in 
the  obtained  values  of  the  wavelength. 

According  to  Geisler's  theory,  a  decay  in  the  direction 
perpendicular  to  the  storm  track,  with  an  e-folding  distance  equal  to  the 
radius  of  deformation  should  be  expected.  In  our  case,  the  value  of  the 
radius  of  deformation  is  less  than  20  kilometers.  On  one  hand,  this 
explains  the  narrowness  of  the  wake  behind  the  storm.  On  the  other  hand, 
as  the  length  scale  of  the  cyclone  is  greater  than  the  length  scale  of 
the  baroclinic  radius  of  deformation,  small-scale  dispersive  wave  modes 
are  not  excited.  Therefore,  a  small  amount  of  energy  is  dispersed 
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Figure  5.  a)  Iscanomal ies  of  the  zonal  velocity  field,  u;  b)  the 

meridional  field,  v;  c)  the  height  field  of  the  top  layer 
for  the  cyclone  travelling  from  west  to  east,  after  7  days 
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towards  the  northern  and  southern  boundaries.  Figures  5a,  b  and  c 
confirm  these  results. 

Permanent  baroclinic  ridges  are  observed  between  upwelled  regions 
(Fig.  5c).  Geisler's  linear  theory  establishes  that  the  amplitude  of 
these  ridges  depends  on  the  vorticity  input  of  the  cyclone  and  not  on  the 
rate  of  input.  The  particular  curvature  of  the  ridges  may  be  explained 
by  the  southward  trajectory  of  the  water  parcels  that  we  have  previously 
discussed. 

The  counterclockwise  motion  of  the  cyclone  generates  a  divergent 
Ekman  flow  towards  the  southern  and  northern  boundaries.  Therefore,  a 
belt  of  "high  water  mass"  is  generated,  by  continuity,  at  both  sides  of 
the  front  (Fig.  5c).  The  water  displaced  away  from  the  front  is  replaced 
by  a  deeper  and  cooler  water  mass,  thus  generating  an  upwelling  effect. 

From  Fig.  5c  it  can  be  seen  that  the  vertical  amplitude  of  the 
upwelling  is  of  the  order  of  10  meters.  Let  us  define  a.H  to  be  such  an 
amplitude  (Fig.  6): 

AH(t)  =  hx(t)  -  hi(o)  (60) 

As  we  are  dealing  with  an  incompressible  fluid,  the  vertically  integrated 
continuity  equation  for  the  top  layer  is: 


1  3  ,  ,  3u  3  v 

- AH  3  (—  +  — ) 

h-  3t  3 x  3y; 


(61) 


where  -  — (<iH)  represents  the  velocity  of  upwelling. 
3 1 


By  simple  scaling,  it  is  inferred  that 

3h 

-fv  3  -g*  — •  +  Tx/p0hi 

3  X 


(52) 
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Figure  6.  Relative  positions  of  the  first  interface,  represented  by 
solid  line,  near  the  front  in  the  (y,  z)  plane  at  times 
t  *  it  and  t  =  0.  The  dotted  line  represents  the  vertical 


motion  of  the  first  interface,  while  the  dashed  I’p.e 


reoresents  the  southern  motion  of  this  same  interface. 
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3h 

fu  s  -g*  ~  +  Ty/p0hx 

3y 

Cross  differentiating  (62)  and  (63)  yields 

,3u  3 v  ,3iy  3tX. 

f(—  +  — )  3  ( - )/Pohl 

3x  3  y  ax  3y  01 

Eliminating  the  divergence  between  (61)  and  (62)  we  have 


(63) 


(64) 


—  (AH)  =  curl  t/p0  f 
3t 


(65) 


where  — (aH)  is  positive,  because  the  curl  of  the  wind  stress  is  positive 
3t 

A  suitable  estimation  of  AH  is  then 


AH  = 


t 

/0  curl  t  dt 


(56) 


By  simple  kinematics,  dt  *  dx/S,  where  S  represents  the  cyclone's  speed. 
Therefore  (66)  may  be  expressed  as 

f  curl  t  dx 


AH  = 


>ofS 


(67) 


By  introducing  the  values  used  in  the  model,  an  approximate  value  of  AH 
is  10  meters.  This  estimate  of  AH  agrees  with  the  values  observed  in 
Fig.  5c. 

The  linear  inviscid  momentum  equations  have  the  form: 


3  u 

—  3  fv  -  g*3h/3x 

at 


(68) 


3  V 

—  *  -fu  -  g*3h/3y 
3t 


(59) 
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The  perturbation  field  is  decomposed  into  a  mean  meridional 
component,  denoted  by  (“)  and  a  wave  meridional  structure,  denoted  by 
(  ).  It  is  observed  from  figures  5a,  b  and  c  that  the  u  and  h  solutions 
are  in  phase  and  that  the  u  and  v  solution  are  in  quadrature.  Thus,  the 
perturbation  field  is  chosen  to  have  the  following  solutions: 

u  *  U(y )  +  U (y )  cos  ij>  (70) 

v  =  Y(y)  +  V(y)  sin  y  (71) 

h  =  TT(y )  +  H(y)  cos  <p  (72) 

where  the  phase,  <|/,  is  equal  to  kx  -  tut.  The  zonal  wavenumber  ,  k,  and 
the  frequency,  w,  are  both  positive  quantities  for  a  perturbation  moving 
eastward  at  the  wake  of  the  cyclone.  Introducing  (70)  -  (72)  into  (68) 
and  (69)  yields 


-<*>  Usirty  -  f(Y+Vsimj>)  =  g*  kHsirnj; 

(73) 

dH  'dH 

in  Vcosip  +  f(U+Ucosip)  =  -g*—  -g*-—  cosv 

dy  dy 

(74) 

Subtracting  the  mean  values,  over  a  period  of  the  wave 

,  allows  us  to 

decompose  (73)  and  (74)  into  the  following  expressions 

• 

7  =  0 

(75) 

— U  -  fV  =  g*kH 

(76) 

r- 

fJJ  =  -g*_ 

dy 

(77) 

dH 

cdV  +  fll  =  -g*—  cos<|> 
dy 

(78) 

As  shown  by  (75),  the  v  perturbation  solution  has  no  mean  meridional 


structure.  This  is  confirmed  in  Figure  (5b)  where  only  the  wave 
meridional  structure  is  observed.  Because  the  minimum  value  of  h  is 
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observed  at  the  front,  its  meridional  derivative  is  zero.  Therefore: 


dH 

dy 


=  0 


(79) 


From  (78)  it  is  inferred  that  V  and  U  have  opposite  signs.  Therefore,  it 
may  be  concluded  that  the  u  perturbation  solution  leads  the  v 
perturbation  solution  by  a  quarter  of  a  period.  On  the  other  hand, 
combining  (78)  with  (76)  yields: 

f2  - 

«(-*  -  1)  !)  *  g*  k  H  (80) 

As  f  is  greater  than  <u,  both  coefficients  of  U  and  H  are  positive.  It 
fallows  immediately  that  the  u  perturbation  solution  is  in  phase  with  the 
h  perturbation  solution.  This  result  is  confirmed  in  Figs.  5a  and  5c. 

It  is  interesting  to  have  an  idea  of  the  order  of  magnitude  of  the 
southward  drift  of  the  subarctic  front.  Fig.  6  gives  us  a  schematic 
representation,  in  the  (y,  z)  plane,  of  the  initial  position  of  the  front 
and  its  location  after  a  time  At.  By  definition: 


,y-W/2, 

3  ^jo  *  Ahjtanh  ; — - — } 

y 


(si) 


Therefore 


dh, 

-  y-W/2  =  — 
<*y  Lv 


dh 


(32) 


This  value  of  —  is  exactly  equal  to  tane.  On  the  other  hand,  from 
dy 

rig.  6  it  is  deduced  that  the  value  of  tane  is  equal  to  the  quotient  of 
the  amplitude  of  the  vertical  displacement  of  the  interface, 


approximately  10  meters,  by  the  southward  shift  of  the  front.  By  simple 
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calculations,  it  follows  that  the  value  of  the  southward  shift  is  of  the 
order  of  10  kms.  In  this  particular  case,  this  is  precisely  the  value  of 
the  width  of  the  front. 

At  the  end  of  the  sixth  day  the  center  of  the  cyclone  is  located  in 
the  middle  of  the  region  of  interest.  Figure  7  shows  the  wave  forced  by 
the  storm's  passage  propagating  towards  the  eastern  boundary.  Tlhe  most 
interesting  feature  to  be  observed  is  that  the  amplitude  of  the  wave  is 
exactly  the  same  at  all  depths.  Furthermore,  no  phase  shift  with  depth 
is  observed. 

Geisler's  theory  predicts  that  very  far  to  the  rear  of  the  storm, 
after  the  dispersion  of  the  inertial  gravity  waves,  the  motion  will  be  in 
geostrophic  balance.  This  is  understandable.  In  the  absence  of 
externally  driven  forces,  the  kinetic  energy  in  the  ocean  will  be 
redistributed.  This  feature  is  observed  nine  days  later.  However,  one 
feature  remains  unalterable:  the  upward  lifting  of  the  four  interfaces 
(Fig.  8).  It  is  expected  that  this  upwelled  region  will  remain  for 
several  weeks  after  the  passage  of  the  storm. 

A  way  of  quantifying  the  dissipation  of  kinetic  energy  after  the 
storm's  passage  is  devised.  For  this  purpose,  the  function 

0  =  D0  exp(-t/x)  (33) 

is  approximated  by  the  least  square  method.  Here  0  represents  the 
maximum  depth  variation  of  the  first  interface;  D0,  the  initial 
amplitude;  t,  time,  and  t,  the  e-folding  time  scale.  The  constants  to  be 


determined  by  the  least  square  method  are  tnD0  and  -l/i.  Values  of  0  are 
evaluated  daily  from  the  seventh  to  the  fifteenth  day.  The  values  of  inD 
versus  time  are  plotted  (Figure  3).  The  resulting  value  of  t  is  of  the 
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Decay  in  time  of  the  amplitude  of  the  vertical  oscillation 
induced  at  the  lee  of  the  cyclone's  track.  While  the  dot 
points  represent  the  original  values  of  In  D  (see  text), 
the  solid  line  represents  the  curve  obtained  by  using  tne 
least  square  method  in  (83). 
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order  of  10  days.  This  relatively  short  e-folding  time  scale  confirms 
the  fact  that  because  of  the  absence  of  any  external  atmospheric  wind 
forcing,  the  kinetic  energy  in  the  ocean  is  rapidly  redistributed.  A 
shorter  e-folding  time  scale  is  to  be  expected  if  vertical  mixing  were  to 
be  included  in  the  model. 

Although  the  lee  side  of  the  storm  is  upwelled,  a  region  of  relative 
maximum  convergence  is  observed  in  the  lower  layers  below  a  region  of 
maximum  divergence  in  the  upper  layer  (Figure  10a,  b,  c  and  d) .  The 
physical  explanation  of  this  phenomenon  is  quite  simple.  The  water  that 
is  upwelled  in  the  top  layer  tends  to  diverge  to  the  meridional 
boundaries.  Therefore,  the  cooler  waters  from  the  deeper  layers  tend  to 
rise  to  the  sea  surface,  generating  a  region  of  relative  maximum 
convergence  in  the  lower  layers. 

The  initial  position  of  the  front  in  the  (y,  z)  plane  is  depicted  in 
Figure  11.  While  the  water  to  the  south  of  the  figure  represents  the 
transition  zone,  the  water  to  the  north  of  tne  figure  represents  the 
subarctic  zone.  The  subarctic  front  represents  the  natural  boundary 
between  these  two  zones.  In  the  transition  zone,  located  at  the  south  of 
the  front,  the  water  is  warmer  and  more  saline  that  in  the  subarctic 
zone. 

In  our  simulation,  the  center  of  the  cyclone  is  travelling  over  the 
subarctic  front  from  west  to  east.  At  the  end  of  the  seventh  day  two 
features  are  observed:  the  southward  drift  of  the  front  and  a  local 
uplifting  of  the  density  interface  (Fig.  12).  Because  of  the  3 -effect, 
the  cross-front  velocity  is  more  intense  at  the  southern  part  of  the 
cyclone.  Thus,  the  Ekman  flow  divergence  will  be  lesser  at  the  northern 
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Figure  10.  Isoanomalies  of  the  height  field  a)  for  the  first;  b)  the 
second;  c)  the  third;  and  d)  for  the  fourth  layers,  for 
the  case  of  the  cyclone  travelling  from  west  to  east 
after  10  davs. 
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Figure  12.  Cross  section  of  the  subarctic  f 
for  the  first  cyclonic  case  afte 
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part  of  the  front.  This  phenomenon  explains  the  slight  southward  shift 
of  the  subarctic  front.  This  particular  form  of  upwelling  is  generated 
by  the  geometrical  form  of  the  cyclone.  An  intensification  of  the  front 
is  also  observed. 

On  the  other  hand,  the  bump  observed  in  Fig.  12  is  generated  by  an 
upwelling  of  different  origin.  This  upwelling  is  driven  by  the  divergent 
Ekman  flow  generated  by  the  counterclockwise  motion  of  the  cyclone. 

Using  (65),  an  attempt  to  estimate  the  height  of  the  bump  is  made. 
Recalling  that  S  is  the  storm  speed,  an  approximate  estimate  of  the 
bump  height,  Ah,  will  be 

PaCDvo2 

Ah  =  - 

Sfp0  (34) 

=  20  meters 

Thus,  the  bump  caused  by  the  Ekman  flow  has  an  approximate  height  of  20 
meters. 

In  general,  a  decrease  of  temperature  is  observed  at  the  wake  of  the 
cyclone.  The  magnitude  of  the  temperature  decrease  is  directly  related 
to  both  the  magnitude  and  the  speed  of  the  cyclone.  3ecause  our  model 
does  not  allow  vertical  mixing,  the  fresher  and  cooler  waters  from 
beneath  cannot  reach  the  sea  surface.  Therefore,  the  temperature  field 
of  the  top  layer  will  remain  almost  unaltered.  However,  a  slight 
southward  shift  of  the  salinity  and  the  temperature  fronts  in  the  upper 
layer  is  detected  nine  days  after  the  cyclone's  passage. 

To  study  the  effects  caused  by  the  latitudinal  change,  a  cyclone 
travelling  from  north  to  south  is  considered.  To  make  comparisons 
easier,  the  stress  and  speed  of  this  cyclone  will  remain  the  same  as  in 
the  previous  case. 
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In  the  present  case,  the  outer  basin  has  an  extension  of  5900  and 
1000  kilometers  in  the  meridional  and  latitudinal  directions, 
respectively.  On  the  other  hand,  the  inner  mesh  has  a  width  of  1000 
kilometers  and  a  length  of  200  kilometers.  The  northern  inner  boundary 
is  located  1500  kilometers  to  the  south  of  the  northern  outer  boundary. 

Initially  located  at  5900  kilometers  north  of  the  equator,  the 
cyclone  travels  along  the  y  axis  at  a  speed  of  500  kilometers  per  day. 

It  reaches  the  northern  inner  boundary  at  the  end  of  the  second  day. 

At  the  end  of  the  fourth  day,  the  center  of  the  cyclone  is  located 
in  the  middle  of  the  southern  inner  boundary.  Essentially,  the  principal 
features  are  recovered,  as  energy  remains  concentrated  in  a  somewhat 
wider  strip  at  the  wake  of  the  storm  (Fig.  13).  This  phenomenon  should 
be  attributed  to  the  increasing  value  of  the  baroclinic  radius  of 
deformation,  which  is  inversely  proportional  to  f. 

It  is  observed  that  the  region  of  maximum  uplifting  is  biased  more 
to  the  right  and  lags  the  center  of  the  cyclone  by  a  distance  of  almost 
300  kilometers.  The  amplitude  of  the  water  uplifted  is  greater,  in  the 
present  case,  by  almost  a  factor  of  2.  This  result  is  consistent  with 
Geisler's  (1970)  theory  in  the  sense  that  it  predicts  a  displacement  of 
the  thermocline,  in  the  wake  of  the  cyclone,  inversely  proportional  to 
f3. 

The  wave  along  the  track  has  a  wavelength  of  500  km,  which  is 
greater  than  in  the  previous  case.  Again,  this  result  is  also  consistent 
with  Geisler's  theory,  which  predicts  that  the  wavelength  should  be 
inversely  proportional  to  f. 
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An  east-west  cross  section  of  the  subarctic  front  one  week  after  the 
cyclone  had  passed  the  front  is  shown  in  Figure  14.  The  main  feature, 
i.e.,  water  uplifted  by  20  meters,  is  observed  in  the  four  layers. 
However,  the  maximum  of  water  uplifted,  in  this  case,  is  slightly  biased 
to  the  west,  i.e.,  to  the  right  of  the  storm's  track.  This  result  is 
perfectly  consistent  with  what  is  observed  in  the  previous  figure.  The 
downwelling  observed  at  both  lateral  boundaries  may  be  explained  by 
continuity.  As  upwelling  persists  one  week  after  the  cyclone's  passage, 
the  result  previously  obtained  still  holds;  namely,  that  in  the  absence 
of  any  external  dissipation,  water  will  remain  upwelled  for  several 
weeks. 

On  the  other  hand,  a  closer  look  at  a  north-south  cross  section  of 
the  subarctic  front,  shows  that  the  features  observed  for  the  cyclone's 
case  previously  considered  are  recovered  (Fig.  15),  namely;  an  upwelling 
of  around  20  meters  of  the  four  interfaces,  a  strengthening  of  the 
subarctic  front,  and  a  bump  of  around  twenty  meters  at  the  site  of  the 
front.  In  the  present  case,  it  is  observed  that  the  southward  drift  of 
the  front  is  more  pronounced  (of  the  order  of  25  kilometers). 

The  oceanic  frontal  response  to  a  severe  atmospheric  front  is  also 
considered.  The  atmospheric  front  is  initially  located  at  1000  km  to  the 
west  of  the  northwest  corner  of  our  region  of  interest  (Fig.  16).  With 
an  angle  of  inclination  of  45°  with  respect  to  the  latitudinal  circles, 
the  atmospheric  front  moves  in  a  direction  perpendicul ar  to  its  axis  of 
orientation.  Travelling  at  a  speed  of  500  kilometers  per  day,  it  reaches 
the  southwest  corner  after  51  hours.  Although  the  northwesterly  winds 
simulated  behind  the  atmospheric  front  are  severe,  no  significant 
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Figure  14.  Cross  section  of  the  subarctic  front  in  the  (x,  z)  plane 
for  the  second  cyclonic  case.  The  solid  line  represents 
the  initial  position  of  tne  density  interfaces,  while 
the  dashed  line  represents  the  position  of  the  interfaces 
after  10  days. 
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15.  Cross  section  of  the  subarctic  front  in  the  (y,  z)  plane 
for  the  second  cyclonic  case.  The  solid  line  represents 
the  initial  position  of  the  density  interfaces,  while  the 
dashed  line  represents  the  position  of  the  interfaces 
after  10  days. 
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Figure  16.  Isoanomalies  of  the  height  field  of  the  top  layer  for  the 
atmospheric  front  case  after  10  days. 
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change  in  the  height  field  is  reported  after  10  days  (Fig.  16). 

Therefore,  it  can  be  concluded  that  the  passage  of  an  atmospheric  front 
plays  no  significant  role  in  the  dynamics  of  the  subarctic  front. 

The  kind  of  atmospheric  wind  forcing  that  is  implemented  in  this 
study  is  the  one  that  is  typically  observed  at  the  latitude  of  the 
subarctic  front.  In  the  first  instance,  the  oceanic  front  is  chosen  to 
be  parallel  to  the  latitudinal  circles.  In  order  to  make  this  study  more 
comprehensive,  an  initial  undulatory  front  is  also  simulated  (Fig.  2). 

At  first  glance  one  would  expect  a  similar  response  of  both  fronts  to 
similar  atmospheric  wind  forcing. 

For  the  first  type  of  oceanic  front,  whenever  the  westerly  wind  is 
simulated,  upwelling  is  observed  at  the  site  of  the  front.  Conversely, 
whenever  a  negative  curl  stress  is  simulated,  convergence  is  observed  at 
the  site  of  the  front.  But  perhaps  the  most  exciting  experiment 
represents  the  simulation  of  a  moving  cyclone.  A  strong  upwelling,  that 
remains  for  several  weeks,  is  observed  along  the  storm  track.  The 
results  of  this  particular  case  agree  very  well  with  Geisler's  linear 
theory. 

The  main  result  of  this  research  is  that  the  momentum  transfer,  in 
the  case  of  the  subarctic  front,  is  confined  to  the  first  150  meters. 

The  only  exception  is  represented  by  the  travelling  cyclone.  Therefore, 
a  good  agreement  with  Roden's  observations  is  achieved. 

The  existence  of  undulatory  fronts  is  very  frequent  in  the  North 
Pacific  Ocean.  Therefore,  an  initial  undulatory  subarctic  front  is 
simulated  by  a  sinusoidal  wave  along  the  latitudinal  circles  (Fig. 

The  longitudinal  wavelength,  Lx,  is  chosen  to  be  500  kms,  while  the 
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meridional  wave  amplitude,  Ly,  is  chosen  to  have  a  value  of  10 
kilometers.  As  our  intention  is  to  observe  the  waves  propagate,  no 
atmospheric  forcing  is  imposed  in  the  first  model  run.  After  ten  days, 
the  undulatory  front  remains  unaltered. 

The  phase  speed  of  a  Rossby-gravity  wave,  cr,  is  given  by 

CR  =  k2+£2+f2/g*H  (35) 

With  the  values  used  in  the  numerical  model,  cr  will  have  an 
approximate  value  of  0.1  cm/sec.  Therefore,  the  westward  propagation 
of  the  Rossby  wave  will  be  of  the  order  of  one  kilometer  in 
ten  days.  This  value  explain  the  fact  that  no  displacement  of  the 
front  is  observed. 

Because  the  meridional  extent  of  the  front  is  small,  the  Rossby  wave 
cannot  propagate  at  an  appreciable  rate  (see  (85)).  It  can  be  concluded 
that  the  3-effect  does  not  play  any  significant  role  in  the  dynamics  of 
the  subarctic  front.  Therefore,  the  f-plane  approximation  made  in 
section  3  is  perfectly  consistent. 

To  make  this  study  more  complete  the  spectrum  of  experiments  is 
broadened.  The  question  to  be  answered  is  whether  a  different  initial 
configuration  of  the  oceanic  front  gives  the  same  kind  of  dynamical 
response  as  in  the  case  of  the  nonundul atory  oceanic  front.  For  this 
purpose,  the  san,c  atmospheric  wind  forcings  described  before  are 
simulated  in  the  case  of  the  undulatory  front. 

The  response  of  the  undulatory  oceanic  front  to  the  passage  of  a 
strong  cyclone  is  addressee.  The  intention  is  to  study  the  oceanic 
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response  of  two  different  initial  configurations,  using  the  same  cyclonic 
wind  forcing.  To  accomplish  this  task,  only  the  height  field  solution  is 
considered  (Fig.  17).  Comparison  between  figures  17  and  10a  shows  that 
the  height  field  variations  are  identical  in  both  fronts.  The  numerical 
simulation  of  the  other  atmospheric  wind  forcings,  i.e.,  westerly  winds, 
curl  and  atmospheric  fronts,  shows  identical  results  in  the  isoanomal 
height  field  for  both  fronts. 

It  can  be  concluded  that  the  frontal  response  to  the  different 
atmospheric  wind  forcings  is  independent  of  the  initial  configuration  of 
an  oceanic  front. 


KILOMETERS 


Isoanomalies  of  the  height  field  of  the  upper  layer  of  the 
undulatory  front  for  the  case  of  the  cyclone  travelling 
from  west  to  east  after  10  cays. 


5.  SUMMARY  AND  CONCLUSIONS 


A  comprehensive  theoretical  and  numerical  study  of  the  subarctic 
front  has  been  performed.  On  the  theoretical  side,  it  has  been 
established  that  horizontal  friction  plays  no  significant  role  as  a 

s 

mechanism  of  dissipation.  A  stability  criterion  has  been  established  by 
which  if  fu-c|2>g*hu‘2/4f2  everywhere,  the  front  is  barotropical ly 
stable.  In  another  theoretical  attempt,  the  numerical  instability  caused 
by  the  constancy  in  the  density  jump  between  layers  has  been  addressed. 

The  numerical  aspect  of  this  research  consisted  of  a  series  of 
numerical  simulations  of  different  atmospheric  wind  forcings  over  the 
subarctic  front.  It  has  been  established  that  the  effects  of  the 
atmospheric  wind  forcing  are  confined  to  the  first  150  meters.  These 
results  agree  witn  Roden's  observations  in  the  sense  that  the  dynamical 
response  to  the  atmospheric  forcing  is  limited  to  the  layer  between  the 
sea  surface  and  the  high  stability  layer  (which  in  the  case  of  the 
subarctic  front  has  an  annual  variation  between  120  and  130  meters  of 
<epth).  The  only  exception  is  represented  by  the  passage  of  a  strong 
cyclone,  where  an  upwelling,  of  the  order  of  20  meters,  of  the  four 
interfaces  is  observed  at  the  wake  of  the  cyclone.  The  following 
conclusions  for  this  case  can  be  drawn: 

1)  The  e-folding  time  scale  of  tne  maximum  depth  variation  of 
the  first  interface  is  of  the  order  of  10  days.  Thus,  in  the 
absence  of  external  forcing,  the  kinetic  energy  is  rapidly 
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2)  The  water  remains  upwelled  for  several  weeks. 

3)  An  intensification  of  the  subarctic  front  takes  place. 

4)  A  southward  shift  of  the  subarctic  front  is  observed. 

5)  A  region  of  relative  maximum  convergence  is  observed  in  the 
Tower  layers  below  a  region  of  maximum  upwelling  in  the  upper 
1 ayer. 

6)  The  u  and  v  fields  are  in  quadrature  in  space  with  the  v  field 
1 eading  the  u  field. 

7)  Energy  is  concentrated  to  a  narrow  strip  at  the  wake  of 
the  storm.  The  result  agrees  with  oeisler's  linear  theory. 

The  simulation  of  the  westerly  atmospheric  winds  shows  a  southward  drift 
of  the  front  due  to  the  southerly  Ekman  flow.  While  the  passage  of  an 
atmospheric  cold  front  does  not  play  any  significant  role  in  the 
dynamical  structure  of  the  subarctic  front,  the  simulation  of  the  curl 
shows  the  classical  convergence  that  strongly  favors  frontogenesis  in 
that  region. 

Two  initial  forms  of  oceanic  fronts  are  considered.  The  results,  in 
all  the  cases  tested,  are  very  conclusive:  the  frontal  response  to  the 
different  atmospheric  wind  forcings  is  independent  of  the  initial  form  of 


the  subarctic  front. 


APPENDICES 


A.  ON  THE  EIGENVALUE  PROBLEM 


During  the  course  of  this  research,  it  was  verified  that  if  the  jump 
in  density  between  layers  remains  constant,  the  numerical  model  becomes 
numerically  unstable.  The  problem  is  traditionally  overcome  by  changing 
the  constancy  in  the  density  jump  between  layers  (McNider  and  O'Brien, 
1973).  However,  a  deeper  understanding  is  required.  To  achieve  this 
goal,  the  phase  speed  of  the  first  baroclinic  mode  is  evaluated. 

To  accomplish  this  task,  a  linearized  version  of  the  non-vi scous  model 
equation  is  sufficient.  Such  a  linearized  version  of  the  momentum 
equations  is 
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where  the  superscript  j  goes  from  1  at  the  top  to  5  at  the  bottom,  and 
again  v  is  the  maximum  value  of  the  pair  (i,  j).  The  linear  form  of  the 
continuity  equation  leads  to 


iii  +  hjW  -  =  o 

at  '‘ax  a  y 

where  H'J  represents  the  mean  depth  of  the  j •  t h  layer.  The 
differentiation  of  (A.l)  with  respect  to  time  yields 
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Introducing  (A. 3)  into  this  last  expression  yields 

32tf3  _  f 3v^  =  q  y  rr£s-pv]  ra2uk  +  a2uki i  (a.5) 

3t2  at  k-l  pj  3x2  3x3y 

This  set  of  equations  is  coupled.  To  decouple  (A.5),  uJ  and  vJ  can 

be  decomposed  into  its  vertical  modes  in  the  following  manner: 

()j  *  l  C.  3  Ult  (A. 6) 

k=l  k 

where  the  value  k  -  5  corresponds  to  the  barotropic  mode  while  k  =  1,2,3 
and  4  represents  the  four  possible  baroclinic  modes.  On  the  other  hand, 
CfcJ  and  u|<  represent  the  eigenvector  and  the  amplitude  of  the  k-th  mode, 
respectively. 

The  elements  of  a  certain  matrix  N|<j  may  be  defined  in  the  form 

;<kJ  -  gHk  (A.7) 

^  J 

Upon  substitution  of  (A. 6)  and  (A.7)  in  (A.5),  a  single  expression  for 
each  mode  is  obtained: 

c.  f  i!k)  =  l  N.  ■  c.  +  liik)  (A. 8) 

*  3t2  3t  k*l  *  3x2  3x3y 

The  set  of  equations  is  no  longer  coupled.  Because  c^J  is  an 
eigenvector  of  the  matrix  N|<j,  it  must  satisfy  the  relation 

Ji"kJ  ck'  ■  V*  (A-9) 

where  a  is  the  associated  eigenvalue.  Using  (A. 9)  in  (A. 8)  a  final 
k 

expression  for  the  linear  u  momentum  equation  is 

3^k  f  3*k  .  „  k  +  j2u,<; 

3t2  at  k  3x2  axsy' 


(A. 10) 


Proceeding  in  the  same  manner,  a  similar  expression  may  be  derived 
for  the  v  momentum  equation: 

^  +  f  auk  =  (3fvk  +  j  (A. 11) 

at2  at  k  axay  ay2 

From  (A. 7)  it  follows  that  the  elements  N5 j  and  Nk5  are  identically 
zero.  This  result  should  be  no  surprise  as  the  barotropic  mode  has  been 
filtered  out.  It  is  convenient  to  define  a  matrix  Mkj  which  is  related 
to  the  matrix  t\j  by  the  form 

Mkj  =  Nk_i  j.i  (A. 12) 

where  k  and  j  vary  from  one  to  four  in  Mkj. 

The  next  task  is  to  evaluate  the  eigenvalues  of  the  first  baroclinic 
mode.  Two  cases  are  considered.  In  one  case  the  densities  are 
arbitrarily  chosen  to  be  pi  =>  1.000,  p2  =  1.002,  P3  »  1.0022,  ?4  =  1.0025 
and  P5  =  1.0027  g  cm"3.  In  the  other  case,  the  density  jump  between 
layers  is  constant  from  layer  to  layer.  Thus,  pj  =  1.000  g  cm"3, 

Pl+i  =  pi  +  iApj  where  i  =  1(1)4  and  Ap^  *  10"3.  In  both  cases,  the 
layer  depths  are 

=  mk 

H2  =  nk  (A. 13) 

h3  =  h4  -  k 

where  k  =  50  meters  and  m  and  n  are  allowed  to  vary  al ternatively  from 
one  to  three.  In  all  the  different  tests  conducted,  the  eigenvalues  of 
the  first  baroclinic  mode  proved  to  be  larger  in  the  case  of  equal  jump 
in  density.  That  difference  ranged  between  30%  and  140%.  (Table  3) 
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The  phase  speed  is  equal  to  the  square  root  of  the  eigenvalue  of  the 
respective  baroclinic  mode.  Therefore,  the  phase  speed  of  the  first 
baroclinic  mode  is  larger  for  the  case  where  the  jump  in  density  is 
constant.  It  is  immediately  conclusive  from  the  CFL  criterion  for 
numerical  instability  that  a  smaller  time  step,  at,  for  this  case  is 
required.  This  result  was  confirmed  in  all  the  model  runs  tested.  In 
particular,  the  time  step  used  in  the  numerical  model  is  20%  larger  than 
the  one  which  would  be  needed  if  an  equal  density  jump  between  layers 
were  used. 
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Table  3 


H 

H.2 

1 

2 

3 

1 

14388.9 

14859.0 

15886.9 

2 

27678.5 

27800.8 

27973.7 

3 

41124.3 

41191.5 

41273.7 

(a) 

H2 

H1 

1 

2 

3 

1 

34534.3 

44855.7 

57612.9 

2 

50452.5 

56025.1 

64769.9 

3 

68917.0 

72054.1 

77099.3 

(b) 

Eigenvalues 

of  the  first  baroclinic 

mode  for  ^  =  mK,  and 

H2  =  nK,  where  m  and 

n  varies  from  1 

to  3  and,  when:  a) 

P1  =  1,0»  p2 

=  1.002 

;  p3  s  1  .0022;  p 

=  1.0025;  p5  =  1.C027 

g  cm-3;  b)  p 

l+i  3  pi 

+  iApj ,  where  i 

=  1(1)4  and  APi  =  10"3 

g  cm" 


B.  ON  THE  BAROTROPIC  INSTABILITY  OF  THE  SUBARCTIC  FRONT 


The  phenomenon  by  which  eddies  are  generated  by  horizontal  shear  in 
the  mean  flow  is  known  as  barotropic  instability.  In  this  case,  eddies 
may  extract  energy  from  the  mean  flow  kinetic  energy.  However,  the  mean 
flow  may  also  extract  energy  from  the  mean  potential  energy  field. 
Therefore,  in  the  case  of  barotropic  instability,  energy  may  be  supplied 
by  the  mean  potential  energy  field. 

Kuo  (1949)  studied  the  stability  of  a  horizontal,  rotating, 
non-divergent  flow.  In  doing  so,  he  extended  Lord  Rayleigh's  (1880) 
theorem  for  a  non-rotating  system.  Kuo  proved  the  necessary  condition 
for  barotropic  instability,  by  which  the  gradient  of  absolute  vorticity 
must  change  sign  at  least  once  along  the  meridional  extension  of  the 
basin  being  considered. 

The  aim  of  the  present  work  is  to  derive  a  stability  criterion  for  a 
geostrophically  balanced  upper  oceanic  front.  To  accomplish  this  thrust, 
a  two-layer  model  is  considered,  in  such  a  way,  that  the  lower  layer 
is  chosen  to  be  infinitely  deep.  To  simplify  the  complexity  of  the 
equations  to  be  used,  the  following  hypotheses  are  made: 

1)  the  oceanic  front  is  geostrophically  balanced;  i.e., 

a77  _ 

-fu  =  g*—  ,  v  =  0,  where  g*  represents  the  reduced  gravity  and 

ay 

the  over-bar  (~)  quantities  represent  the  basic  state; 

2)  because  the  meridional  extent  of  the  oceanic  front,  Ly,  is  of 
the  order  of  thirty  kilometers,  the  r-plane  approximation  is 
used; 
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3)  the  zonal  wavenumber,  k,  is  much  smaller  than  (Ly)_1;  therefore 
only  long  waves  are  considered; 

4)  the  vorticity  at  the  front,  represented  by  -7u/3y,  is  much  less 
than  the  planetary  vorticity,  f. 

The  linear  inviscid  perturbation  equations  for  a  single-layer  fluid 
may  be  written  as 
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where  the  prime  quantities  represent  the  perturbation.  In  (3-2)  the 
geostrophic  balance  is  substracted. 

Assuming  that  the  perturbation  quantities  have  a  wavelike  form 


u' 

| 

u(y) 

v' 

* 

v(y) 

exp  [ik(x-ct)] 

(B-4) 

h' 

H(y) 

where  U(y),  V(y)  and  H(y)  are  respectively  the  amplitudes  of  u',  v1  and 
h‘ ,  the  above  set  of  equations  has  the  form 

du 

ik(u-c)U  +  (  —  -  f J V  =  -g*  ikH  (3-5) 

dy 

fl!  +  ik(U-c)  V  =  -g*  dH/dy  (B-6) 

iknU  +  —  (TTV)  +  ikCu-c)  -  0  (3-7) 

dy 

This  system  of  eauations  may  be  reduced,  after  some  algebraic 
manipulations,  to  the  single  equation 
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fh'  (hu* ) ' 

(hH* )  *  -  [  —  +  — -  +  a2T?]  H  =  0  (B-8) 

U-C  U-C 

where  a2  =  f2 ( 9*^) “1  is  the  inverse  of  the  Rossby  radius  of  deformation 
squared. 

Since  the  meridional  length  scale  we  are  interested  in  is  of  the 
order  of  thirty  kilometers,  it  may  be  assumed  that  H(»)  =  H(-~)  =0. 

For  this  type  of  perturbation  to  be  unstable,  the  phase  speed, 
c(=  cr+  i c -j ) ,  must  be  complex.  The  amplitude  functions  must  also  be 
compl ex. 

For  the  sake  of  simplicity,  the  variable  Z(=  TT-c)  is  introduced,  and 

l 

the  following  change  of  variables  is  made  (Howard,  1961):  H  =  ZT  Q  . 

With  these  considerations,  (8-8)  may  be  expressed  as 

■  |  ^ 

CfiZQ')'  -  [fP  +  1  (TT'TP  ) '  +  a2hZ  +  ~-]Q  =  0  (B-9) 

Multiplying  (B-9)  by  Q*  and  upon  integration  by  parts  yields 

/{PIQ'I2  +  [fP  +  l-  (Pu')'  +  a2hZ  +  I Q | 2 } dy  =  0  (B-10) 

where  the  limits  are  omitted. 

The  imaginary  part  of  Z  is  -c-j.  It  follows  immediately  that  the 
imaginary  part  of  Z_1  is  -Ci|ZJ"2.  Therefore,  the  imaginary  part  of 
(3-10)  yields 

hu'2 

Ci/CKIQ1!2  +  {a2lT  -  — t}  |Q|2]  dy  =  0  (3-11) 

| 

Whenever 

r*hu 1 - 

|Z|2  =  |  u-c | 2 >R2  =  Z-—  (3-12) 
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the  integrand  of  (3-11)  has  a  positive  definite  form.  In  such  a  case, 
(8-11)  holds  only  if  ci  =  0.  If  the  condition  (B-12)  is  satisfied,  then 
for  a  given  phase  speed  c  for  every  point  along  the  front,  the  stability 
of  the  front  to  the  fluctuation  is  guaranteed. 

Looking  at  the  problem  from  a  mathematical  point  of  view,  criterion 
(B-12)  is  met  outside  the  circle  of  radius  R  centered  at  (U,  o)  along  the 
cr  axis.  For  a  given  oceanic  front,  both  "h  and  IT  are  functions  of  the 
latitude,  y.  Therefore,  for  every  latitudinal  point,  y,  across  the  front 
there  corresponds  one  interval  and  one  radius. 

The  two  extremes  -  maximum  and  minimum  -  of  this  family  of  intervals 
delimit  the  region  of  stability.  In  other  words,  any  point  outside  that 
region  corresponds  to  a  stable  fluctuation  of  the  front.  This  is  a 
sufficient  condition  of  stability. 

What  follows  is  an  attempt  to  apply  this  stability  criterion  to  an 
upper  oceanic  front.  For  the  sake  of  generality,  let  us  define  the  depth 
of  the  upper  layer,  "FT,  via  a  hyperbolic  tangent,  in  such  a  way  that 

"FT  =  "FT00  -  Ah  tanh(y/Ly)  (3-13) 

where  "FT00(Ah)  represents  the  mean  value  (amplitude)  of  the  upper  layer. 
For  this  steady  oceanic  front  to  subsist,  the  interface,  IT,  must  be 
balanced  initially  by  a  zonal  velocity  of  the  form 

U  =  -(g*/f)(~J  (3-14) 

By  continuity,  the  meridional  velocity,  v,  is  initially  set  to  be  zero. 

It  is  convenient  here  to  non-dimensional ize  the  governing  equations. 
Let  y  be  scaled  by  Ly,  w  by  f,  k  by  a  and  U  and  c  by  the  zonal  velocity 


at  the  front,  noted  |ujm.  The  non-dimensional  variables  n,  <,  x>  Y  and 
are  introduced  as  follows: 


n  =  y/Ly  k  =  u/f 

Y  =  c|U|m  x  =  k/a 
w  ■  "u/Rm 

The  basic  state  may  be  rewritten  as 

"R  =  h0n  -  Ah  tanh  n 


—  9*  Ah 

u  -  —  —  sech-  n 

'  ~y 


(B-15) 


(B-16) 


( B— 17) 


The  expression  for  the  zonal  velocity  at  the  front  has  the  form 


(B— 18) 


By  definition,  it  follows  that  u  will  be  equal  to  secl^n.  Thus, 
the  derivative  of  y  with  respect  to  n  has  the  form  -2  tanh  n  secifn. 
With  this  formulation,  criterion  (3-12)  has  the  final  form 


9*^00  Ah 

[ u -Y 1 2 >-t- -y  ( 1  -  - —  tanh  n)  sech4  n  tanh2  n  =  $2 
f^Lu^  ri 

y  oo 


(3-19) 


For  the  purpose  of  clarity,  a  family  of  circles  of  radius  each 
one  centered  at  u,  is  constructed  in  the  complex  y  plane  (Fig.  13).  The 
envelope  defined  by  the  family  of  circles  delimits  the  outer  region  of 
stability.  Because  of  the  positive  contribution  of  the  term  hZ | Q ' | 2  in 
(B-ll),  it  would  be  inappropriate  to  state  that  the  region  delimited  by 
the  points  inside  the  envelope  corresponds  to  the  unstable  case. 

Let  yQ  and  y,  be  the  minimum  and  maximum  real  values  of  the  points 
along  the  envelope,  respectively.  Three  cases  are  possible: 
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Circumference  of  equation  jyr  -  a1-4- 
comolex  v  plane  for  :nj«l. 
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(B-20a) 

(B-20b) 

(B-20c) 


(i)  0<y0<y1 

(11) 

(Hi)  y0<0<y1 

(i)  If  both  Yg  and  yl  are  positive,  the  stability  in  the  x  complex  plane 
corresponds  to  the  points  outside  the  interval  (x;  *  x0)  since  the 
dimensionless  frequency  '<  is  positive  and  real  (Fig.  19a).  ( i  1 )  The  same 

result  holds  if  both  extremes,  yq  and  yl»  are  negative  (Fig.  19b). 

(iii)  If  y0  is  negative  and  y1  is  positive,  the  stability  region  is 
determined  by  the  points  inside  the  Interval  (xQ ,  xx)  (Pig*  19c). 

For  h00  =  100  m.  Ah  =  10  m,  Ly  =  10  km,  y0  and  y1  are  found  to  be  0 
and  1.25  m  sec-1,  respectively.  For  a  two-day  period,  these  values 
correspond  to  a  wavelength  less  than  21.6  kms.  Thus  a  sufficient 
condition  of  stability  for  a  typical  upper  oceanic  front  is  established. 
This  criterion  guarantees  that  waves  whose  wavelength  are  larger  than  22 
kms  may  be  barotropically  stable. 


C.  LIST  OF  SYMBOLS 


a  coefficient  of  thermal  expansion 
A  horizontal  diffusion  coefficient  of  momentum 
b  coefficient  of  saline  contraction 
B  positive  constant  defined  by  Eq.  (15b) 
c1,  c2  arbitrary  constants,  defined  by  Eq.  (54) 
cr,  Ci  real  and  imaginary  parts  of  the  phase  speed 
c  phase  speed 

Cg  surface  drag  coefficient 
D  depth  variation  of  the  first  interface 
Do  amplitude  of  the  depth  variation  of  the  first  interface 
E  Ekman  number 
f  Coriolis  parameter 
F(y)  left  hand  side  of  Eq.  (45a) 
g  gravity  acceleration 
g*  reduced  gravity 
G(y)  left  hand  side  of  Eq.  (45b) 
h  thickness 
H  mean  height 
k:  zonal  wavenumber 

!<  positive  constant,  defined  by  Eq.  (26) 

1:  meridional  wavenumber 

Lx,  Ly  latitudinal  and  meridional  frontal  length  scales 
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Mkj»  Nkj  matrices  of  elements  {k,  j)  defined  by  Eqs.  (A-7)  and  (A-12), 
respectively 
p  pressure 

Pa  atmospheric  pressure 
Q  function  defined  by  Eq.  (8-9) 
r  cyclone's  radii 

r0,  vo  positive  constants,  defined  by  Eq.  (23) 

R  positive  constant,  defined  by  Eq.  (3-12) 

Ro  thermal  Rossby  number 
S  velocity  of  translation  of  the  cyclone 
5  salinity 
t  time 

T  time  scale  of  dissipation 
7  temperature 

T0,  Tl  positive  constant  amplitudes  of  the  wind  stress,  defined  by 
Equations  (22)  and  (25),  respectively 
u  zonal  velocity  component 
v  merdional  velocity  component 
V  meridional  mass  transport 

Vc,  Vp  general  and  particular  solutions  of  the  meridional  mass 
transport,  V 
W  width  of  the  basin 
X  constant  defined  by  Eq.  (2H) 

z  independent  vertical  coordinate,  increasing  upward 
Z  function  defined  by  Eq.  (3-9) 
a  inverse  of  the  Rossby  radius  of  deformation 
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8  Coriolis  parameter 

Y i  n,  u,  <»  x  non-dimensional  variables,  defined  by  the  system  (B-15) 

5  free  surface  of  the  uppermost  layer 
ah  amplitude  of  the  interfaces'  depth  variation 
ah  bump  height 
AH  upwelling  amplitude 
A21  temperature  amplitude 

e  angle  formed  by  the  upwelling  and  the  southward  drift  of 
the  oceanic  front 

0  angle  defined  by  Equations  (22a  and  22b)  respectively 
\  zonal  wavelength 
u  maximum  value  of  the  pair  (i,  j) 
pa  air  density 
Pj,  p0  sea  density 
a  eigenvalue 

tx,  horizontal  wind  stress  components 

u  circumferential  velocity  of  the  cyclone 
T  e-folding  time  scale 
x,  c  arbitrary  functions  defined  by  Eq.  (55) 
iii  phase 

'f  forcing  function  defined  by  Eq.  (57) 
u  frequency 
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